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Section  1 


BACKGROUND  AND  OBJECTIVES 

LIGHTNING  INDUCED  VOLTAGE  IN  AIRCRAFT  ELECTRICAL  CIRCUITS 


One  of  the  trends  in  the  design  of  modern  aircarft  is  toward  the  use  of  minia- 
turized solid-state  electronics  in  avionics,  automatic  flight  control,  and  other  func- 
tions. The  decrease  in  weight  and  power  consumption  which  these  devices  afford  has 
enabled  improvements  in  performance  and  economy,  so  the  trend  is  likely  to 
become  widespread  in  the  design  of  new  military  and  civilian  aircraft.  Because 
of  the  inherently  small  size  and  low  operating  power  levels  required  by 
miniaturized  solid-state  electronics,  however,  these  components  have 
been  found  to  be  more  vulnerable  to  electromagnetic  interference  than  their 
vacuum  tube  counterparts  of  an  earlier  day.  Such  interference  may  result 
from  on-board  systems  such  as  radio  transmitters  or  relay  operation,  or 
from  external  sources  such  as  lightning  or  nuclear  electromagnetic  pulse 
(NEMP). 

Incidents  have  been  reported,  for  example,  in  which  solid-state  elec- 
tronics have  been  upset  or  permanently  damaged  as  a result  of  lightning 
strikes  to  aircraft  (Refs.  1 and  2).  A number  of  research  programs  have 
been  conducted  by  various  Government  agencies  or  research  laboratories 
to  determine  the  extent  of  lightning  related  interference  and  the  mechanisms 
by  which  it  occurs  in  aircraft  electrical  circuits  (Refs.  3-0),  and  it  has  been 
learned  that  electromagnetic  fields  caused  by  lightning  may  appear  inside 
typical  aircraft  and  induce  transient  surge  voltage  in  electrical  wires  and 
cables.  These  lightning- induced  voltages  are  in  addition  to  those  which  may 
enter  aircraft  electrical  circuits  as  a result  of  direct  lightning  stroke  con- 
tact with  externally  mounted  electrical  components  such  as  navigation  lights 
or  antennas. 

Malfunction  of  sensitive  electronics  may  occur  if  the  induced  voltage 
exceeds  its  overvoltage  withstanding  capability  or  if  the  accompanying  in- 
duced current  surge  results  in  the  dissipation  of  excessive  power  in  semicon- 
ductor junctions.  Since  lightning  electromagnetic  fields  usually  permeate  the 
entire  aircraft,  redundant  systems  are  also  susceptible  and  may  not  provide 
their  intended  backup  capability. 


Heretofore,  most  lightning  protection  design  has  been  for  control  of  the 
direct  effects  of  lightning,  such  as  fuel  ignition  and  structural  strength 
degradation,  or  directly  conducted  surge  voltages  and  currents  arising  from 
strokes  to  navigation  lights  or  antennas.  However,  the  increasing  depen- 
dence of  critical  navigation  and  flight  control  functions  on  solid-state  micro- 
electronics has  resulted  in  recognition  of  the  need  for  protection  against 
the  indirect  effects  of  lightning. 
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Test  techniques  and  equipment  have  been  designed  for  subjecting  com- 
plete aircraft  to  simulated  lightning  strokes  so  that  the  degree  of  suscepti- 
bility of  various  individual  circuits /systems  can  be  determined.  Such  tech- 
niques allow  measurement  of  actual  induced  voltages  in  these  circuits  for 
comparison  with  known  component  withstand  levels,  to  determine  the  need 
ior  additional  protective  measures.  To  avoid  expensive  retrofit  programs, 
however,  lightning  protection  should  be  designed  into  each  aircraft  system 
from  the  start.  This  means  that  designers  must  have  information  concerning 
the  expected  susceptibility  of  particular  circuits  while  they  are  still  on  the 
drawing  board,  when  there  is  as  yet  no  aircraft  on  which  to  run  tests. 

To  fulfill  this  need,  a program  was  initiated  by  the  Air  Force  Flight 
Dynamics  Laboratory  at  Wright- Patterson  Air  Force  Base  to  develop  a com- 
puterized analysis  technique  for  calculating  voltages  expected  to  be  induced 
in  typical  aircraft  circuits  by  lightning  stroke  currents  flowing  through  the 
aircraft.  The  overall  objective  was  to  develop  a computer  program,  readily 
and  economically  usable  by  aircraft  designers,  to  assess  the  impact  of  var- 
ious structural  and  electrical  system  design  configurations  on  lightning  sus- 
ceptibility, and  thereby  to  provide  a tool  by  which  design  optimization  and 
tradeoff  studies  can  be  made  from  the  standpoint  of  lightning  protection. 

The  analytical  approach  followed  was  based  upon  a preliminary  attempt, 
made  under  an  earlier  National  Aeronautics  and  Space  Administration  con- 
tract (Ref.  5),  which  showed  promise  when  compared  with  actual  experimen- 
tal measurements.  After  some  introductory  discussions  of  lightning  induced 
voltage  mechanisms,  this  report  describes  the  analytical  steps  applied  and 
the  computer  programs  developed  to  calculate  voltages  induced  in  electrical 
conductors  at  various  locations  inside  a complete  aircraft. 

DIRECT  COUPLED  VOLTAGES 


Directly  coupled  voltages  occur  as  a result  of  direct  contact  of  lightning 
strike  currents  with  exposed  (external)  electrical  assemblies,  such  as  an- 
tennas and  navigation  lights.  If  a lightning  flash  punctures  a lamp  globe  or 
antenna  housing  so  that  direct  contact  may  be  made  with  a filament  or  antenna 
element,  a portion  of  the  lightning  current  may  be  conducted  into  associcted 
electrical  wiring.  This  voltage  will  be  accompanied  by  a voltage  surge 
limited  in  amplitude  by  the  insulation  breakdown  voltage  level  of  the  electri- 
cal assembly  or  associated  wiring,  whichever  is  less.  Unless  external  pro- 
tection is  applied  to  prevent  puncture  of  the  external  assembly,  it  is  often 
damaged  beyond  operational  capability.  Even  if  this  damage  is  acceptable, 
however,  the  surge  voltages  and  currents  which  proceed  into  associated 
wiring  are  usually  hazardous  to  connected  equipment  such  as  power  control 
or  communication  electronics. 

Thus,  protection  against  these  surges  must  be  provided.  The  magnitude 
of  these  conducted  surges  and  the.  adequacy  of  protective  devices  designed  to 
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control  them  must  usually  be  evaluated  by  full-scale  simulated  lightning  tests 
of  the  external  assemblies  in  question.  Government  specifications  for  some 
of  these  devices  or  protective  equipment  are  now  in  existence  (Ref.  7).  Di- 
rectly coupled  voltage  and  current  surges  are  considered  direct  effects 
were  not  dealt  with  further  in  this  program. 

INDIRECT  COUPLED  VOLTAGES 


The  other  mechanism  by  which  lightning  can  affect  aircraft  electrical 
and  avionics  systems  is  by  the  generation  of  magnetically  induced  and  resis- 
tive voltage  rising  within  aircraft  electrical  circuitry.  To  describe  the  man- 
ner in  which  induced  voltages  occur  it  is  first  necessary  to  consider  the 
mechanisms  by  which  magnetic  fields  appear  inside  sin  aircraft. 

For  a long  conductor  carrying  a current,  i,  and  whose  return  path  is 
far  removed,  the  average  field  intensity  at  a distance,  r,  from  the  conductor 
is 


as  shown  in  Figure  1. 


(a)  Current  Carrying  Filament 


H = o ( r < r0) 


(b)  Tubular  Conductor 


(c)  Irregular  Conductor 

Figure  1.  Magnetic  Fields  Around  Current  Carrying  Conductors 
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If  instead  of  a solid  wire  the  current  were  carried  on  a hollow  tube  of 
radius  rQ,  as  shown  on  Figure  1(b),  the  field  intensity  at  radius  r < rQ,  would 
be 


H 


(2) 


2irr 

and  at  the  surface  of  the  tube  where  r equals  rQ  the  field  intensity  would  be 

i 


H = 2„r„ 

Since  the  circumference  of  a tube  is 

C = 2irr0 


(3) 


(4) 


it  follows  that  the  field  intensity  at  the  surface  of  a tube  is 


(5) 


The  average  current  density  at  the  surface  of  the  tube  is  also  equal 
total  current  divided  by  the  circumference: 

_ i 

^aw  ~n 


(6) 


If  the  conductor  is  not  cylindrical,  as  shown  in  Figure  1(c),  the  field  intensity 
at  different  points  on  the  surface  will  be  different.  Field  intensity  will  still 
be  equal  to  the  total  current  divided  by  the  circumference: 


The  actual  field  intensity  will  be  greater  than  average  at  points  where  the  ra- 
dius of  curvature  is  less  than  average,  and  less  than  average  at  points  where 
the  radius  of  curvature  is  greater  than  average,  as  shown  in  Figure  2. 


For  example,  in  a wing  carrying  lightning  current,  the  leading  and 
trailing  edges  have  radii  of  curvature  much  smaller  than  average.  Field  in- 
tensity along  the  leading  and  trailing  edges  should  then  be  quite  high  compared 
to  the  field  intensity  along  the  top  or  bottom  surfaces. 

Since  both  the  average  current  density,  JAvE,  and  the  average  field  inten- 
sity, HavE,  are  equal  to  the  total  current  divided  by  the  circumference, 

I - H -- 

°AVE  nAVC  C 

(8) 

it  follows  that  the  tangential  field  intensity  at  the  surface  of  a conducting  ob- 
ject is  equal  to  the  current  density  at  that  point.  This  is  in  fact  true,  at 
least  for  transient  currents.  The  relation  is  not  true  for  d-c  currents  or 
transients  sufficiently  slow  that  appreciable  magnetic  fields  penetrate  the 
skin. 
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H < H,v, 


Figure  2.  Field  Intensity  Versus  Radius  of  Curvature 

The  orientation  of  the  H field  vector  is  always  at  right  angles  to  the  di- 
rection of  the  current  vector,  as  shown  in  Figure  1. 

While  small  gaps  in  the  structure  direct  the  current  around  the  gap,  the 
magnetic  field  is  virtually  unaffected,  except  directly  on  the  surface  and  on  a 
length  scale  that  is  small  compared  to  dimensions  of  the  gap  interrupting  the 
current  flow  (see  Figure  3). 

So  far,  only  the  field  external  to  the  aircraft  has  been  dealt  with. 
Even  if  the  aircraft  has  an  electrically  continuous  metallic  skin,  some  mag- 
netic flux  can  appear  within  the  aircraft  as  lightning  current  diffuses  through 
the  skin  to  the  inside  surface.  Cancellation  effects  will  eliminate  this  flux 
in  perfectly  symmetrical  cases,  such  as  a cylinder  with  uniform  skin  current 
distribution,  but  in  other  cases  some  net  diffusion  flux  may  exist  inside.  The 
interior  field  is  generally  characterized  as  having  a slower  time  to  crest  than 
the  exterior  field  as  well  as  a lower  amplitude;  this  is  illustrated  in  Figure  4. 

If  apertures  exist  in  the  aircraft  skin,  a portion  of  the  external  magnetic 
flux  surrounding  all  of  the  current  flowing  through  the  aircraft  will  leak  in- 
side through  these  apertures,  as  shown  in  Figure  5.  This  is  known  as  aper- 
ture  flux;  and  it  appears  inside  much  sooner  than  the  diffusion  flux,  since  its 
velocity  is  unimpeded,  and  has  a higher  rate  of  change,  sinilar  to  that  of 
the  total  lightning  current.  Aperture  flux  is  usually  more  localized  than  dif- 
fusion flux,  and  in  the  vicinity  of  apertures  it  may  have  a much  higher  ampli- 
tude than  the  diffusion  flux. 
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Aileron 


Figure  3.  Current  Flow  and  Magnetic  Field  Around  Structural  Gaps 


H„t  (Step  Pulee) 


Figure  4.  Internal  Diffusion  Fields  as  a Function  of  External  Fields 
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Figure  5.  Aperture  Fields 

An  aperture  can  be  described  in  terms  of  an  equivalent  magnetic  dipole 
(Figure  6).  In  Figure  6(a),  current  streamlines  are  seen  being  diverted 
around  a hole  in  a current  carrying  sheet.  A filamentary  dipole  producing 
the  same  magnetic  effects  as  the  diverted  current  flow  would  be  as  shown 
on  Figure  6(b).  The  magnetic  field  pattern  produced  by  such  a dipole  is  the 
same  as  the  classic  magnetic  field  produced  in  the  near  field  zone  by  a mag- 
netic dipole,  and  is  shown  in  Figure  7.  The  farther  one  is  from  the  opening, 
the  less  is  the  field  intensity,  decreasing  approximately  as  the  third  power  of 
the  distance,  for  distances  that  are  large  compared  to  the  size  of  the  opening. 

The  changing  internal  magnetic  fields  link  electrical  wires  and  cables 
inside  the  aircraft,  inducing  voltages  therein.  The  induced  voltages  are  re- 
lated to  the  lightning  current  by  inductive  transfer  functions  (Ref.  3)  in  accor- 
dance with  Faraday's  law.  Since  the  airframe  is  composed  of  inactive  cir- 
cuit constants,  the  transfer  function  for  diffusion  flux  coupling  should  be  a 
constant  inductance  for  any  lightning  waveform,  relating  the  portion  of 
lightning  current  appealing  at  the  inside  surface  of  the  skin  to  the  voltage  it 
induces  in  a circuit.  The  transfer  function  relating  voltages  induced  into  a 
circuit  by  the  aperture  flux  is  probably  more  complex  because  of  less  uniform 
field  patterns  and  aperture  geometries. 

In  addition  to  the  magnetically  induced  voltages,  the  resistance  of  the 
metallic  skin  will  permit  resistive  voltage  differences  in  the  skin  (or  struc- 
ture) along  the  path  of  lightning  current  flow.  If  an  aircraft  electrical  cir- 
cuit employs  the  structure  as  its  return  path,  then  this  resistive  voltage 
enters  the  circuit,  in  series  with  the  magnetically  induced  voltage  in  the 
same  circuit  and  any  other  (normal)  steady- state  operating  voltages  present. 
Capacitively  coupled  voltages  may  also  be  produced  in  these  circuits;  however, 
the  essentially  uniform  conducting  skin  of  metallic  aircraft  keeps  potential 
differences  among  structural  elements  low,  thereby  limiting  the  voltages 
which  can  be  electrostatically  coupled  to  interior  electrical  circuits.  In 


(a)  Actual  Current  Flow- 


Current  Flow 
in  Equivalent 
Dipole  Antenna 


Current 

Streamlines 


(b)  Equivalent  Dipole 


Figure  6.  Development  of  the  Equivalent  Magnetic  Dipole 


Field  Intensity 
Varies  with  Distance 
from  Dipole 


Figure  7,  Field  Pattern  Due  to  an  Aperture  Dipole 
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practice,  experimental  measurements  have  shown  magnetic  and  resistive 
components  to  be  the  most  predominant  (Refs.  3-6). 


The  combination  of  the  resistive  and  magnetic  components  of  induced 
voltages  should  therefore  be  expressible  as  follows: 

dr(l-e_oct)  yt)]  + M diL  (9) 


eoc  = R.i^t)  + 


dt 


a~ST 


where: 

eoc  = voltage  induced  in  the  circuit  (using  the  airframe 
as  return) 

Rs  = effective  structural  resistance 

Mj  = diffusion  transfer  inductance  between  lightning  current 
flowing  on  the  inside  surface  of  the  skin  and  the  particu- 
lar electrical  circuit 

Ma  = aperture  transfer  inductance  between  the  total  lightning 

current  flowing  through  the  aircraft  and  the  particular  aix'- 
craft  circuit 

iL (t) - lightning  current  (a  time  varying  function) 

a = reciprocal  of  the  time  constant  of  current  penetration  into 
the  aircraft  skin 

Of  course,  circuit  transmission  line  and  termination  impedance  characteris- 
tics as  well  as  secondary  induced  effects  may  change  the  induced  voltage 
appearing  at  a particular  point  from  that  predicted  by  Equation  9.  Equation  9 
is  therefore  most  appropriately  viewed  as  representing  the  induced  source 
voltage  driving  the  distributed  aircraft  circuit. 

in  the  first  experimental  programs  (Refs.  3-5),  transfer  functions  de- 
rived from  induced  voltage  data  indicated  that  most  of  the  enclosed  magnetic 
flux  was  of  diffusion  origin,  and  the  M?term  of  Equation  9 was  not  necessary 
for  this  equation  to  adequately  describe  the  measured  induced  voltage  wave- 
forms. The  work  of  Refs.  3,  4,  and  5 was  conducted  on  an  E89.1  fighter 
aircraft,  however,  which  has  few  apertures.  Subsequent  work  on  different 
aircraft  (lief.  6 ) with  more  apertures  showed  evidence  of  much  greater 
aperture  field  coupling  into  aircraft  circuits;  this  mode  was  often  more  pre- 
dominant than  either  the  diffusion  magnetic  or  resistive  mechanisms.  At 
the  conclusion  of  the  F89J  tests,  work  was  initiated  on  a completely  analyti- 
cal technique  to  arrive  at  the  same  transfer  functions  (Ref.  5).  This  involved 
a mathematical  representation  of  an  E89J  wing  and  an  electrical  circuit  con- 
ductor inside.  Some  simplifying  assumptions  relating  to  wing  geometry  and 
lightning  current  flow  were  made  in  this  attempt,  and  the  magnetic  flux  linking 
the  conductor  and  its  airframe  was  calculated  as  a function  of  an  assumed 
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lightning  current  filament  in  the  wing  skin.  The  contributions  from  a large 
number  of  such  filaments,  assumed  to  comprise  the  wing,  were  summed  to 
obtain  the  total  magnetic  flux  linking  the  conductor  and  its  airframe  return. 
From  this,  the  transfer  inductance,  M1#  was  derived.  The  resistive  trans- 
fer function,  Rs,  was  calculated  as  a function  of  geometry  and  material  resis- 
tivity. The  resulting  values  of  Rs  and  M,^  compared  well  with  corresponding 
values  derived  from  measured  induced  voltages  on  a circuit  inside  the  F89J 
wing.  The  work  accomplished  in  this  program,  particularly  that  dealing 
with  diffusion  coupled  voltages,  is  based  on  this  preliminary  approach. 

PROGRAM  OBJECTIVES 


The  basic  objective  of  this  program  was  to  develop  computerized  analyt- 
ical models  to  determine  possible  lightning  induced  voltages  in  aircraft  elec- 
trical circuits.  Specific  requirements  for  these  models  were  that  they  repre- 
sent the  major  airframe  sections  of  a complete  aircraft,  including  fuselage, 
wing,  horizontal  stabilizer,  and  vertical  stabilizer. 

Another  goal  was  to  incorporate  as  many  refinements  over  the  original 
model  of  Reference  5 as  possible.  The  desired  improvenents  included: 

1.  Calculation  of  the  actual  lightning  current  distribution  throughout 
the  circumference  of  each  major  section.  (The  original  model 
assumed  a uniformly  distributed  current.  ) 

2.  Representation  of  circuit  conductors  of  different  lengths  than  that 
of  the  major  airframe  section  itself. 

3.  Location  of  the  circuit  conductor  anywhere  inside  the  airframe, 
instead  of  along  its  axis  of  symmetry  only. 

4.  Calculation  of  voltages  induced  by  aperture  flux,  such  as  would 
penetrate  holes,  windows,  and  access  doors.  (The  original 
model  assumed  a completely  enclosed  airframe.  ) 

5.  Calculation  of  the  effect  of  varying  one  circuit  location  or  airframe 
geometry  parameter  while  holding  the  other  constant. 

To  the  extent  possible  within  the  program  resources,  it  was  also  desired 
to  represent  internal  structural  elements,  such  as  spars,  ribs,  and  bulkheads, 
and  to  accommodate  more  complex  electric  circuit  configurations  such  as 
shielded  cables,  wire-to-wire  (independent  return),  and  individual  circuit 
impedance  characteristics. 

Upon  completion  of  each  basic  model,  its  mathematical  equations  and 
validity  were  to  be  verified  by  having  them  represent  simple  geometries  for 
which  textbook  solutions  are  available  and  aircraft  geometries  for  which  test 
data  are  available. 
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The  computerized  models  were  to  be  programmed  in  FORTRAN  extended 
version  IV  for  execution  with  punched  cards  on  the  U.  S.  Air  Force  CDC  6600 
computer  at  Wright- Patterson  Air  Force  Base,  Ohio.  The  program  was  to 
be  delivered  as  a punched  card  deck.  A user's  manual  and  a final  technical 
report  were  also  to  be  delivered. 

BASC  APPROACHES 


As  previously  discussed,  lightning  induced  voltages  in  aircraft  electrical 
circuits  occur  because  time  varying  aperture  and  diffusion  magnetic  fluxes 
exist  inside  the  airframe.  Aperture  flux  penetrates  through  openings  such  as 
windows  and  access  doors  in  the  aircraft  structure.  Diffusion  flux  appears 
inside  as  lightning  currents  diffuse  through  the  thickness  of  the  metallic  skin 
and  appear  on  its  inside  surface.  • 

Because  the  methods  by  which  these  fluxes  enter  the  airframe  are  fun- 
damentally different,  it  was  decided  that  completely  separate  models 
should  be  developed  to  represent  the  diffusion  and  aperture  coupling  mechanisms. 
The  diffusion  model  is  based  on  the  original  approach  of  Reference  5,  which 
assumed  no  apertures,  whereas  the  aperture  model  is  based  on  treatment  of  a 
single  aperture  which  opens  into  a relatively  small,  confined  space  in  the  airframe. 
Contributions  to  aperture  flux  from  individual  apertures  are  considered  of 
greatest  interest,  because  the  flux  entering  from  one  aperture  is  frequently 
segregated  from  that  entering  through  other  apertures  by  the  presence  of 
spars,  ribs,  bulkheads,  and  ? .her  interior  walls,  which  act  as  electromagnetic 
shields. 

For  most  airframe  or  circuit  situations  it  is  not  intuitively  obvious  which 
of  these  two  fluxes  induces  the  greater  voltages.  Therefore,  it  will  be  neces- 
sary for  designers  to  utilize  both  models  for  a complete  evaluation  of  possible 
induced  voltages;  but  as  experience  is  gained,  situations  will  become  apparent 
which  heavily  favor  one  or  the  other  model. 
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Section  2 


DIFFUSION  MODEL 


INDUCE*  VOLTAGE  THEORY 


For  the  diffusion  model  it  was  necessary  to  relate  lightning  currents  flowing 
in  the  aircraft  skin  to  voltages  induced  in  aircraft  electrical  circuits  inside . 

To  derive  this  relationship,  two  fundamental  laws  were  utilized: 

• The  Biot-Savart  law,  which  describes  the  density  of  magnetic  flux  at 
a specific  point  away  from  a current  carrying  conductor. 

• Faraday's  law,  which  describes  the  voltage  induced  in  a conductor  by 
a changing  magnetic  flux  passing  through  a loop  formed  by  this  con- 
ductor. 

Any  loop  formed  by  an  electrical  conductor  such  as  an  aircraft  electrical 
circuit,  which  is  linked  by  a changing  magnetic  field,  will  have  voltage  induced 
in  it  equal  to  the  negative  time  rate  of  change  of  the  total  magnetic  flux  linking 
the  loop.  This  is  Faraday's  law  and  is  expressed  as: 

e =-£!  (10) 

m dt 

where: 


em  = total  emf  (volts) 
dr  = total  flux  (webers) 
t = time  (seconds) 

It  was  next  necessary  to  relate  the  total  flux,  Y,  to  the  lightning  current. 

The  magnetic  flux  which  links  an  open  surface  such  as  that  surrounded 
by  an  aircraft  electrical  circuit  (including  its  return  path)  can  be  found  by 
integrating  the  flux  density,  B,  over  the  surface  area  linked  by  B.  This  may 
be  expressed  as 


* - JJ,  B-ds 


(11) 


where : 


i|f  = total  flux  (webers) 

B = flux  density  (Wb/m2) 
o 

S = surface  area  (m  ) 
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Equations  10  and  11  relate  induced  voltage  flux  to  total  flux,  and  total 
flux  to  flux  density.  In  relating  flux  density  to  lightning  current  it.  is  appro- 
priate to  consider  the  physical  situation  which  exists  when  lightning  strikes 
an  aircraft  structure.  Shown  in  Figure  8 is  a filamentary  representation  of 
an  aircraft  win®,  inside  which  is  located  an  electrical  conductor  and  its  air- 
frame return,  forming  the  circuit  loop  ABCD.  Because  of  the  short  time  dura- 
tion of  most  lightning  strokes,  nearly  all  of  the  lightning  current  flows  in  the 
skin  rather  than  through  internal  spars  and  ribs;  therefore  only  the  skin  is  rep- 
resented. 


Figure  8.  Circuit  Wire  in  Aircraft  Structure 
That  Has  Been  Struck  by  Lightning 

Assuming  that  this  is  so  and  that  lightning  current  flows  in  a lineal  di- 
rection, the  aircraft  structure  can  be  represented  by  a large  number,  n,  of 
parallel  skin  current  filaments.  The  voltage,  V,_D,  appearing  at  the  inboard 
end  and  the  outboard  end  of  the  loop  is  equal  to  the  line  integral  of  voltage  in- 
duced around  t.’ve  loop  ABCD.  As  previously  discussed,  the  voltage  induced 
in  the  loop  is  dependent  upon  the  magnetic  flux  passing  through  this  loop. 

This  vlux  is  in  turn  a function  of  flux  density,  as  indicated  by  Equation  11. 

For  the  arx  angement  shown  in  Figure  9,  the  magnetic  flux  density  produced 


r 


Figure  9.  Current  Carrying  Filament 

at  some  point,  p,  with  respect  to  a current  filament  is  defined  by  the  Biot- 
Savart  law  as 
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where: 


b -■&-  r ^ ai 

n 4?r  J r3 


In  = current  (amperes) 

Bn  = magnetic  flux  density  (Wt/m2 ) 

l,r:  dimensions  (meters) 

p = peimeability  of  the  medium  (for  air  = 47 r x 10  H/m) 

Note  that  at  this  point  an  expression  for  the  flux  density  has  been  intro- 
duced which  is  dependent  upon  current.  Since  the  structure  has  been  represen- 
ted by  a parallel  array  of  n current  carrying  filaments,  there  will  be  n contri- 
butions to  the  flux  density  at  point  p and  all  other  such  points  in  space. 

It  still  remains  to  express  the  flux  density  B in  terms  of  the  airframe 
geometry.  Figure  10  shows  a typical  skin  current  filament  and  the  aircraft 
circuit  loop  previously  considered.  To  obtain  the  total  flux  passing  through 
the  loop  it  is  necessary  to  integrate  the  flux  density  over  the  loop  area.  Equa- 
tion 12  is  expressed  in  terms  of  the  geometry  of  Figure  10.  From  this  is 
obtained 


— r*  r . 1 

JU-zP-i r~  (f-zP+r3  dz  1st  integral  (13) 


+ r ‘ ± 

47r  J ■V{z-i)3‘trs  (z-i)3+r2 


2nd  integral  (14) 
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Figure  10.  Skin  Current  Filament  and  Flux  Density  Through  Aircraft 
Electrical  Circuit  Loop  ABCD 
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The  first  integral  (Equation  13)  can  be  rewritten  as  Equation  15  and  integrated 
by  basic  integral  no.  173  (Ref.  8,  p.  71)  as  follows: 


1st  integral  = ^ J0  (jg8+  z*  . 2S.Z  + rs 

_jA_  [" 2rjj 

4ir  [4ia- 4(ra+ ia)] 

1st  integral  = ~ . ===== L=rl 

[ r J(Z  - i)a+  raJ , 


_jA_  2r(2Z  - 21) 

4jr  [4ia-  4(ra+  ia)]  y Za-  21Z 


+ ra+  iaJ, 


where,  for  the  basic  integral  no.  173  (Ref.  8), 

a = 1,  b = -2i  and  c = (ra+  jfa) 

The  second  integral  (Equation  14)  can  be  rewritten  and  solved  in  the  same 
manner  as  the  first  integral: 


2nd  integral  . <£  ^ 


2iZ  + r-)a 


(Z  - i) 
r /(Z  - i)4 


-] 


It  is  seen  that  the  integral  of  Equation  17  is  evaluated  in  the  Z direction  from 
the  bottom  of  the  filament  at  0 to  i,  and  the  integral  of  Equation  19  is  evalu- 
ated from  S.  to  the  top  of  the  filament  at  L.  This  integration  gives  B as  a 
function  of  position  in  terms  of  £ and  r: 


B«.r)=£ 


i 

r fFTT1 


(L  - i) 
r J(L  - jO5 


Equation  20  is  thus  a general  expression  for  the  flux  density,  B,  at  a point 
at  some  distance  from  any  of  the  current  carrying  filaments. 

Now  that  an  analytical  expression  has  been  developed  for  flux  density, 
the  flux  linking  the  circuit  loop  can  be  determined  by  integrating  the  flux 
density  in  the  manner  suggested  by  Figure  11  and  Equation  21: 

* = J7S  B*  ds  (21) 

The  circuit  loop  inside  an  aircraft  structure  need  not  run  the  entire 
length,  L,  of  that  structure.  It  may  instead  begin  at  any  arbitrary  point,  il. 
and  terminate  at  any  arbitrary  point,  ia.  Therefore  is  the  lower  limit  and 
£ais  the  upper  limit  of  the  first  integration  over  the  circuit  loop  length.  The 
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Figure  11.  Skin  Current  Filament  and  Adiacent  Circuit  Loop 


flux  linking  the  loop  is  also  dependent  upon  the  distance  to  the  loop  location 
from  the  current  filament.  Thus  the  flux  linking  the  loop  in  the  radial  direc- 
tion is  simply  all  of  the  flux  out  to  a radial  distance  D2  minus  all  of  the  flux 
out  to  a radial  distance  D1.  Equation  21  now  has  limits  and  can  be  expressed 
as 


(22) 


Inserting  the  expression  for  B and  performing  the  double  integration 
yields  the  flux  linking  the  circuit  loop  due  to  a single  skin  current  filament. 
This  expression  is  presente<’  as  Equations  23  through  26: 


» = BaL 

4 7T 


/iFT5-!,  io&(S±ZLi*) 


yii 


(23) 

(24) 


( \Ui+  ra  - £1  logs  rJ  + 

( JUi-  L)a+  ra  - (i1  - L)  logc(^  jrljj  (25) 

-)) 


L)a  + ra  - (f  2-  L)  logi 


|V(lg-  L)a+  ra 


+ (j?2  - L) 


11, 


(26) 


D, 


A cross-sectional  view  of  the  situation  shown  in  Figure  11  might  appear 
as  shown  in  Figure  12.  From  this  figure  it  is  clear  that  the  circuit  loop  need 
not  be  in  the  same  plane  as  the  skin  current  filament.  Note  that  the  time 
varying  current  that  forms  a part  of  Equations  23  through  26  makes  tlv>  flux 
a time  varying  function,  as  required  by  Faraday's  law  (Equation  10). 


The  flux  linking  the  circuit  loop  caused  by  the  current  filament  can  be 
calculated  by  assigning  appropriate  values  to  and  t2,  computing  the  value 
of  V when  r = l)2,  and  then  subtracting  from  this  the  value  of  * when  r = Dj. 
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Figure  12.  Cross-Sectional  View  of  Wing  Showing  Distances  Used 
to  Compute  Magnet  Flux  Passing  Through  Circuit  Loop 
ABCD  by  Equations  23-26 

The  total  flux,  which  links  the  loop  due  to  all  current  filaments,  is  the 
summation  of  fluxes  from  all  such  filaments  that  pass  through  the  same  loop. 
Since  all  filaments  will  be  at  different  distances  Dj  and  D2  from  the  loop,  the 
evaluation  of  Equations  23  through  26  must  be  performed  n times.  The  transfer 
inductance,  M,  between  a conductor  carrying  a current  and  another  circuit 
through  which  flux  generated  by  the  first  conductorpasses  is  generally  defined  as 

M -tioui  (27) 

* Tot «] 

The  total  transfer  inductance  is  therefore  the  sum  of  all  of  the  fluxes 
for  all  filaments,  divided  by  the  total  current  responsible  for  that  flux,  or 

n = n 

M = £ y2-  (28) 

n = 1 ItoUI 

This  inductive  transfer  function,  when  inserted  into  Equation  9,  enables  ex- 
pression of  the  magnetically  induced  voltage  in  an  aircraft  electrical  circuit 
as  a function  of  the  lightning  current. 

SKIN  CURRENT  DISTRBUTION  THEORY 

Experimental  measurements  of  skin  currents  in  aircraft  (Ref.  6)  have 
indicated  that  lightning  currents  do  not,  in  fact,  distribute  evenly  around  the 
circumference  of  an  airframe  cross  section.  In  all  but  uniformly  symmetrical 
bodies  (e.  g. , a cylinder)  the  current  in  each  filament  comprising  the  body 
will  be  somewhat  different  from  the  current  in  its  neighbors.  Accordingly, 
a subroutine  was  developed  to  calculate  the  amount  of  current  flowing  in  each 
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of  the  skin  current  filaments  comprising  the  airframe  sections  (Ref.  9).  This  sub- 
routine, which  is  based  on  inductive  current  division,  is  described  in  the  fol- 
lowing paragraph. 

CURRENT  DIVISION 

Figure  13  shows  mutually  coupled  inductances  through  which  current  flows 
and  voltage  is  developed.  If  there  are  two  circuits  (Figure  13a),  then 


Only  the  bilateral  case,  in  which  Mia  = Mai,  will  be  treated  here.  This  is  no 
real  restriction,  because  in  all  physically  realizable  systems  mutual  induc- 
tance is  bilateral.  Only  the  case  in  which  all  currents  are  in  phase  --  the 
low-frequency  case  --  will  be  treated.  While  in  physical  systems  this  need 
not  be  so,  there  are  many  systems  in  which  current  division  is  controlled 
only  by  inductive  effects.  Purely  for  ease  of  numerical  analysis,  only  the 
frequency  for  which  (d/dt)  is  numerically  equal  to  unity  will  be  considered. 

The  analysis  is  valid  for  other  frequencies,  subject  only  to  the  above  restrictions. 


Figure  13.  Mutually  Coupled  Inductances:  a)  Two  Circuits,  b)  N Circuits 


19 


Under  the  above  conditions, 


V i = Ljij  - Miaia 

(31) 

Vs  = -M2iii  + Ljia 

(32) 

In  the  general  case,  Figure  13(b), 

Vj  ~ L ji  j ~ Ml2i2  - M jjij  .... 

“M  in 

(33) 

V2  - "Majij  + Lai2  - M^ia.... 

-M2b  iB 

(34) 

V3  = -M3li1  - M38i2  ■ L3i3  . . . . 
• 

-M3,  ia 

(35) 

• 

vn  = -Mniii  - Mn2ia  - MB3i3  . . . 

• Ejin 

(36) 

Equations  33  through  36  may  be  placed  in  matrix  notation  as 


vx 

Lj  * Mj2  ■ M y , , , , 

-Mia 

ii 

v2 

-M  21+  La  - M 23  . . . . 

-M» 

la 

v3 

• 

= 

-Mai  - M32  - L3  .... 
• 

-Ms, 

X 

ia 

• 

V. 

• 

-Mni  - Mb8  “ Mb3  . . . , 

“Lnn 

• 

in 

or,  in  more  compact  notation: 

I V | --  | M I x I i I 

Multiplying  by  the  inverse  of  the  M matrix,  |M  |_1: 

| M |"1  x | V | = | M I"1  x | M | x | i | 
or: 

| i | = | M I"1  x 1 V | 


(37) 


(38) 

(39) 

(40) 


ii 

i2 

ia 

• 

• 

li 

• 

in 

where  mu,  mu,  mia  . . . . 


mu  m 22  i3  • • • • jn 


m*i  m 22  my  ....  m^ 
m3l  m32  maa  ....  m^ 


mtj  ninj  mB3  • • • • niBB 


x 


Vi 

V2 

V3 

i 


(41) 


are  the  elements  of  the  inverse  of  the  M matrix. 


If  all  of  the  voltages  are  the  same  and  equal  to  V,  as  in  the  case  if  all  of 
the  inductances  are  connected  in  parallel,  the  absolute  current  in  each  element 
is 

h = <mu+  mi2+  mi3 + min)V 

i2  = (m8l+  ma+  m23 + m^)  V 
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(42) 

(43) 


ia  = (m3i  + m38  + m^  . . . . + m3n)V 


(44) 

(45) 


iB  = (mBl  + mB8+  rrins  . . . . + mu)V 

The  total  current  that  flows,  which  is  proportional  to  the  impressed  voltage, 
is 


ii  = (ii  "*■  i3  (3  "*■  • • • • in)V 

The  fraction  of  the  total  current  that  flows  in  each  circuit  is 


3 it 

T -22 

ia  - T 


I =h 

in  • 


(46) 

(47) 

(48) 

(49) 

(50) 


SELF  AND  MUTUAL  INDUCTANCES 


This  analysis  treats  the  case  in  which  the  self  and  mutual  inductances 
are  those  of  parallel  circular  conductors  of  a sufficient  length,  compared  to 
the  spacing  between  conductors,  that  all  end  effects  may  be  ignored.  Only 
the  case  in  which  all  conductors  are  far  removed  from  any  conducting  sur- 
faces such  as  ground  will  be  considered. 


Figure  14  shows  a single  conductor  in  space,  carrying  a current,  I.  The 
magnetic  field  intensity  in  the  space  around  this  conductor  is 


H 


The  magnetic  flux  density  is 

B = /uH 


4?r  x 10"71 
2 tt 


= 2 x 10~7-  Wb/ms 
r 


(51) 


(52) 


The  self- inductance  of  the  conductor,  i,  is  defined  as 

Li=-f^-Wb/A  (53) 

The  ratio  of  webers  per  ampere  is,  of  course,  given  the  name  ’’henries." 

The  magnetic  flux,  cp,  is  equal  to  the  area  under  the  B curve  (Figure  14)  from 
rx  (the  conductor  surface)  out  to  some  other  point,  R,  which  defines  the  return 
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Figure  14.  Self-Inductance 


path  for  the  current,  I,  in  the  conductor.  If  one  postulates  a conductor  car- 
rying direct  current  and  located  in  free  space,  this  return  path  will  be  at 
infinity.  The  flux  density  goes  to  zero  as  R goes  to  infinity,  but  the  area 
under  the  curve,  cp,  also  goes  to  infinity.  If  cp  goes  to  infinity,  then  L,  as 
defined  in  Equation  53,  also  goes  to  infinity.  Accordingly,  one  cannot  speak 
of  a single  value  as  describing  the  inductance  of  an  isolated  conductor. 

If  the  conductor  is  carrying  a transient  or  alternating  current  rather  than 
a direct  current,  an  inductance  can  be  defined;  this  is  because  the  magnetic 
fields  cannot  instantaneously  fill  the  entire  region  around  the  conductor  but, 
instead,  propagate  outward  at  the  speed  of  light.  Because  the  effective  dis- 
tance to  which  they  propagate  is  time  or  frequency  dependent  the  inductance 
will  also  be  time  or  frequency  dependent.  In  this  analysis,  R is  taken  as  the 
distance  to  which  a field  could  propagate  in  one  microsecond  --  300  meters. 
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The  area,  tp,  under  the  B curve  of  Figure  14  is 


<*>  = 2 x io“^  i rR~ 

Jr1  r 

(54) 

tp  = 2 x 10“^  I log€  r J 

R 

Tl 

(55) 

cp  = 2 x 10~7I  log.f 

(56) 

Remembering  the  definition  of  L (Equation  53), 

L,  = f = 2x  10”7  loge 

_R_ 

ri 

(57) 

The  mutual  inductance  between  conductors  i and  j is  defined  as 

(58) 

cp  j.  the  flux  linking  conductor  j and  set  up  by  the  current  I,  in  conductor 
i (as  shown  in  Figure  15),  is 


<Pj  = 2 x 

(59) 

cpj  = 2 x 

10-7  It  loge  ~ 

1 9 

(60) 

I 


Hence, 


MtJ  = 2 x KTloge  f- 


(61) 


PROGRAM  OPERATION 


With  the  cross  section  of  the  airframe  in  the  X-Y  plane,  the  X,  Y,Z  coor- 
dinates of  the  conductors,  and  their  radii,  are  read  and  stored  in  a matrix, 
printed  for  inspection,  and  reread.  The  arbitrary  distance  to  which  the  fields 
propagate  (300  meters)  is  given  as  R5  in  the  computer  listing.  The  spacing 
between  all  conductors  is  then  calculated,  and  the  mutual  inductances  are  cal- 
culated and  loaded  in  the  array.  Self-inductances  are  loaded  into  the  appro- 
priate elements,  those  on  the  main  diagonal. 


At  this  stage  the  matrix  holds  the  absolute  currents  in  the  individual  con- 
ductors (assuming  V = 0),  currents  corresponding  to  those  given  in  Equations 
42  through  45.  The  total  current  is  then  calculated;  the  fractional  current  is 
then  calculated  and  stored. 


COMPUTER  PROGRAM  DIFFUSION 


GENERAL  DESCRIPTION 


The  computer  program  DIFFUSION  was  established  to  represent  an  air- 
craft as  a combination  of  several  independent  sections.  Each  of  these  sections 
is  represented  in  the  computer  program  by  an  array  of  parallel  current  car- 
rying filaments.  Figure  16  shows  a complete  aircraft,  while  Figures  17 
through  20  show  the  individual  sections  of  the  aircraft  modeled  by  this  program. 


Figure  16. 


Complete  Aircraft  Represented  by  Parallel  Current 
Carrying  Filaments 
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Each  section  of  the  aircraft  is  completely  described  by  several  geometric 
constants,  from  which  the  computer  program  calculates  the  location  of  each 
current  filament  with  respect  to  a coordinate  system.  Individual  sections 
are  shown  in  Figures  17  through  20. 

The  geometrical  dimensions  X 1,  Rl,  Yl,  and  Y2,  etc.  are  read  into  the 
computer  program  as  the  first  step  in  execution.  At  the  same  time,  the  initial 
location  of  an  enclosed  electrical  circuit  conductor  and  a set  of  modifiers  are 
read  in.  These  modifiers  allow  the  user  to  reposition  the  electrical  conductor 
during  program  execution. 


The  variations  which  are  made  under  program  control  are  enumerated  and 
illustrated  below,  using  a fuselage  section  as  an  example: 

(l)  The  X coordinate  (Figure  21)  may  be  varied  horizontally  in  a step- 
wise manner  from  X f (initial  X coordinate)  to  Xf  (final  X coordinate). 

0 The  Y coordinate  may  be  varied  vertically  in  a stepwise  manner  from 
Yj  (initial  Y coordinate)  to  Yt  (final  Y coordinate)  (Figure  22). 

(3)  The  length  or  position  of  the  circuit  conductor  may  be  varied  hori- 
zontally in  a stepwise  manner  by  varying  the  Z coordinate  of  either 
or  both  of  the  conductor  end  points,  Zx  and  Zs  (Figure  23). 

0 Any  combination  of  the  X,  Y,  and  Z coordinate  variations  may  also 
be  made.  The  X coordinate  may  be  varied  until  it  reaches  a partic- 
ular value  (0  - 0 ),  after  which  the  Y coordinate  may  be  varied 
( (?)  - 0 ) until  it  reaches  a particular  value;  then  the  Z1  and/or 
Z2  coordinates  may  be  varied  until  a final  position/length  is  achieved 
( ® - ® ) (Figure  24). 

Incrementing  of  all  three  coordinates  may  occur  sequentially,  simulta- 
neously, or  in  combination.  Thus  variation  of  one  variable  need  not  be  ter- 
minated prior  to  changing  the  value  of  another  of  the  variables  (see  Figure  25). 

For  each  circuit  conductor  location  the  program  then  determines  the 
magnetic  flux  density  at  the  forward  or  inboard  end  of  the  circuit  conductor 
as  shown  in  Figure  26.  It  then  computes  the  transfer  inductance  between  the 
circuit  formed  by  the  enclosed  conductor  and  airframe  return  and  the  current 
filaments  used  to  represent  the  aircraft  section  under  investigation. 


Once  the  transfer  inductance  and  resistance  values  have  been  computed, 
the  open  circuit  voltage  versus  time  is  tabulated  for  Equation  62: 


eoe=  Rri(t)+  M 


d(l  - e~a')  i(t) 
dt 


(62) 


The  user  may  input  as  many  different  sets  of  data  cards  as  are  desired. 
The  program  will  execute  each  set  over  the  range  of  values  given  and  print 
out  the  data  generated  for  each  set. 
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Top  View  Fuselage 
Change  of  Posit,  in  OnJ> 


Change  of  Position  and  1 ength 


Top  View 
Change  of  l ength 

Permissible  Variation  of  Z Coordinates  of  Enclosed 
Elect rical  Conductor 


Figure  24.  Permissible  Variation  of  Enclosed  Electrical 
Conductor  Coordinates 


AX  / 0-* 
AY  / Of 
A Z • 0 


AX  / 0 
AY  0 
A 7 0 


Figure  25.  Example  of  a Possible  Set  of  Variations  of  Circuit  Conductor 
Location  and  Length 


f igure  26,  Location  of  Enclosed  Circuit  Loop 
and  Flux  Density  Computation 


TTie  computer  program  initially  divides  the  continuous  geometrical  struc- 
ture into  an  array  of  parallel  current  carrying  filaments  and  computes  the  cur- 
rent distribution  in  each  filament,  using  the  method  already  described  under 
"Skin  Current  Distribution  Theory. " It  then  defines  a horizontal  plane  defined 
by  the  circuit  conductor  and  a return  conductor  in  the  airframe  skin.  The 
program  computes  the  flux  density,  B,  at  the  forward  end  of  the  circuit  con- 
ductor (Figure  26)  and  the  flux  passing  through  the  defined  plane  contributed 
by  each  of  the  current  filaments. 

These  flux  linkages  are  summed  and  divided  by  the  total  lightning  current, 
to  obtain  the  transfer  inductance,  M.  If  the  location  of  the  circuit  conductor 
is  to  be  repositioned,  for  design  optimization  studies,  the  computer  program 
input  data  establish  the  step  size  and  direction  in  which  to  move  the  electrical 
circuit  conductor  for  the  second  operation.  In  such  a case,  new  coordinates 
of  another  horizontal  plane  are  determined  and  the  flux  density  and  flux  compu- 
tations are  performed  again.  Each  time  the  operation  is  performed,  a flux 
density  is  determined  at  the  new  location  of  the  forward  end  of  the  enclosed 
electrical  conductor,  as  well  as  the  total  flux  linking  the  newly  defined  circuit. 
After  each  set  of  conditions  has  been  calculated  the  program  determines  whether 
there  are  other  geometries  to  be  evaluated. 

DIFFUSION  FLOW  DIAGRAMS 

An  elementary  flow  diagram  of  the  DIFFUSION  computer  program  is  shown 
here  as  Figure  27;  Figure  28  is  a detailed  flow  diagram  of  the  program.  A 
program  listing  for  DIFFUSION  is  given  in  Figure  29*;  the  listing  includes,  in 
addition  to  the  main  portion,  subroutines  MATRIX,  MATINV,  and  MATZER. 

The  program  begins  (lines  1-103)  with  some  introductory  comments  to  assist 
the  user  in  operation.  It  next  establishes  five  files  in  which  to  temporarily 
store  data  generated  by  the  program. 

In  lines  170-210  the  program  initializes  the  constants  to  be  used  in  the 
computation  and  reads  in  a control  variable,  A,  which  routes  the  program 
to  line  240,  1840,  1860,  or  1880,  depending  on  the  geometry  specified  by  the 
user.  Upon  selection  of  the  appropriate  geometry,  the  computer  program 
reads  in  the  data  pertinent  to  that  geometry,  and  then  prints  out  a heading 
to  indicate  that  diffusion  coupling  is  being  computed  in  the  particular  geometry 
named.  The  initial  data  read  contains  the  location  of  the  circuit  conductor  to 
be  evaluated  and  a set  of  modifiers  with  which  the  user  may  change  the  position 
of  the  electrical  circuit  conductor  inside  the  particular  geometry.  The  user 
may  make  as  many  modifications  in  these  data  as  he  desires;  for  each  modi- 
fication, one  program  execution  occurs. 

After  the  circuit  conductor  location  has  been  defined,  the  geometry  that 
has  been  selected  can  be  described  as  an  array  of  current  carrying  filaments 
whose  locations  are  computed  from  the  constants  of  geometry  and  the  math- 
ematical expressions  derived  to  analytically  define  that  geometry.  This  is 
done  in  lines  270-1371  for  a fuselage,  or  lines  187-2580  for  wing  type  geom- 

*This  listing  is  for  a program  that  will  be  run  on  the  General  Electric  Time 
sharing  computer.  A program  listing  for  the  CDC6600  machine  is  included 
in  Appendix  III,  "Program  Listings  for  CDC6600  Computer." 
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Figure  27.  Elementary  Flow  Diagram  of  Diffusion  Program 
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Figure  28.  Detailed  Flow  Diagram  of  Diffusion  Program  (Sheet  1 of  23) 
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Figure  28.  Detailed  Flow  Diagram  (Sheet  2 of  23) 
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Figure  28.  Detailed  Flow  Diagram  (Sheet  5 of  23) 
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I Figure  28.  Detailed  Flow  Diagram  (Sheet  6 of  23) 


38 


Writ*  into  HARPO  the 

valttftS 

X » XJ;  Y • Y2 


Compute  coordinates  of  current  fila- 
ments in  upper  left  comer  movinf 
count  ere  lock vise. 

Set  M - 0 


Lines 

1291-1360 


Define  angle  8 - Tl: 
T1  « MD/R1 
X2  * H2-R1  sin  (Tl) 
Y2  * H2 *R  1 cos  (Tl) 


Define  the  nearest  horizontal  plane  which  includes 
the  circuit  conductor  and  structure. 


Yes 


Is  the  conductor 
in  a corner  ? 


Figure  28.  Detailed  Flow  Diagram  (Sheet  7 of  23) 
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Figure  28.  Detailed  Flow  Diagram  (Sheet  10  of  23) 
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Figure  28,  Detailed  Flow  Diagram  (Sheet  16  of  23) 
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Figure  28.  Det 


Figure  28.  Detailed  Flow  Diagram  (Sheet  18  of  23) 
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Figure  28.  Detailed  Flow  Diagram  (Sheet  20  of  23) 
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Figure  28.  Detailed  Flow  Diagram  (Sheet  22  of  23) 
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1C  DIFFUSION A COMPUTER  PROGRAM  WHICH  CALCULATES  THE 

PC  DIFFUSION  FIELDS  AND  THE  DIFFUSION  COUPLED 

3C  VOLTAGES  INTERIOR  TO  SEVERAL  AIRCRAFT 

4C  GEOMETRICAL  COMPONENTS* 

SC 

AC 

1C  KEITH  J.  MAXWELL  BLDG  9-809  GENERAL  ELECTRIC  COMPANY 
8C  100  W0OOLAWN  AVE*  PI TTSFI ELD* MASS*  01801 
9C  PHONE  (41 3) -494-3531* 

IOC 

11C 


I8C  DEVELOPED  UNDER  CONTRACT  F336I I-74-C-3060  USAF  FLIGHT 
13C  DYNAMICS  LABORATORY. 

I4C 

ISC 

I6C  THE  PROGRAM  READS  DATA  FROM  AN  EXTERNAL  FILE  THE  NAME 
I7C  OF  WHICH  HAS  BEEN  SET  TO  "MAXWELL*** THE  INPUT  DATA  SHOULD 
18C  BE  ARRANGED  AS  FOLLOWS  FOR  FUSELAGE  GEOMETRIES! 

I9C 

80 C LINE  NUMBER  100  A 

8IC  110  AI*RI*R8*XI  *Y1*  C7»X5»  D3* D4*  D5*  Y5*  D6*  D7*  DO 

88C  180  S*L7*0*L3*  03*  0 4*  0 5*  L 4*  06*  07*  08*  TO#  T9*  1 4* 

83C  Gl*G8*G3*G4 

84C  130  A 

8SC  140  - SAME  AS  ABOVE  USING  8ND  DATA  SET 

86C  LINE  NUMBERS  MAY  BE  ADDED  INDEFINITELY  UNTIL  ALL  CASES  HAVE 
87C  BEEN  DESCRIBED  * 

SBC 

89C 

30C  DATA  ARRANGEMENT  FOR  WING.H0R1Z  STAB* AND  VERT  STAB 
SIC  SHOULD  BE  AS  FOLLOWS! 

36C 

33C  LINE  NUMBER  100  A 

34C  110  A8» R* Cl*  T*  S*L7*  C7*XS*  03*  D4*  DS*  Y5*  D6*  D7»  D8 

35C  180  O*L3*O3*04*O5*L4«06*07*08*  TB*  T9*  1 4* 

36C  61  *G8»  63*64 

37C  1 30  A 

SBC  140 SAME  AS  ABOVE  USING  8ND  DATA  SET— — 

39C  ADDITIONAL  LINES  OF  DATA  MAY  BE  USED  UNTIL  ALL  CASES  ARE 
ADC  DESCRIBED*  GEOMETRIES  MAY  BE  MIXED  OR  SEPARATED  AS  DESIRED* 
41C 
48C 

43C  A DESCRIPTION  OF  THE  VARIABLES  FOLLOWS! 


44C 

A.  ... 

... . Tur 

17. k nr  M . DaiiTr*  nr  dbad.m  ti  tur 

46C 

GEOMETRICAL  C9NFI GURATION* 

AfC 

Diwri 

^Dv 

H*l  * 

’-■--rusoLHot 

50  C 

Me* 

A«3* 

’••••WIWU 

HORIZ  STAB 
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51 C A*4 VERT  STAB 

52C  THE  VALUES  AI.A2.A3.A4  ARE  USED  AS  A C0MPARI S0N  WITH  THE 
53C  VALUE  OF  A T0  INSURE  THAT  THE  INPUT  DATA  C0RRESP0NDS  T0  THE 
54C  GEOMETRY  SPECIFIED* 

55C  THE  VALUES  0F  R|  AND  R2  ARE  THE  RADIUS  0F  CURVATURE  0F 
56C  THE  T0P  C0RNERS  AND  THE  B0TT0M  C0RNERS  0F  THE  FUSELAGE 
57C  RESPECTIVELY. 

SBC  XI  AND  Y 1 ARE  THE  HEIGHT  AND  WIDTH  0F  THE  FUSELAGE* 

59C  C7  IS  USED  T0  DECIDE  H0W  MANY  RELOCATIONS  OF  A CIRCUIT 
60C  CONDUCTOR  ARE  TO  BE  MADE. 

61 C X5  AND  Y5  ARE  THE  INITIAL  X-Y  COORDINATES  OF  A CIRCUIT 
62C  CONDUCTOR.  THE  CIRCUIT  BEGINS  AT  A DEPTH  OF  L3  INSIDE 
63C  THE  FUSELAGE  AND  EXTENDS  TO  THE  DISTANCE  L4. 

64C  A SET  OF  MODIFIERS  IS  PROVIDED  FOR  EACH  VALUE  DESCRIBING 
6SC  THE  LOCATION  OF  THE  CIRCUIT*  THESE  MODIFIERS  CHANGE  THE 
66C  ORIGINAL  POSITION  OF  THE  CIRCUIT  BY  A STEP  SIZE  GIVEN 

X STEPPED  BY  AN  AMOUNT 

Y -STEPPED  BY  AN  AMOUNT 


67C  AS  t 
6BC 
69  C 
TOC 


L3* 
L 4- 


D3 
D6 

•STEPPED  BY  AN  AMOUNT  03 
•STEPPED  BY  AN  AMOUNT  06 


71 C STEPPING  BEGINS  AT 

72C 

73C 

74C 


E*D4 
E*D7 
E*0  4 

E*07 

75C  FOR  THE  VARIABLES  X.Y.L3.L4  RESPECTIVELY 
76C  STEPPING  OF  ANY  ONE  VARIABLE  TERMINATES  WHEN 


77C 

78C 

79C 

80C 


E*D5* 

E*D8* 

£E«05- 

E»08* 


-X-XMAX 

*Y*YMAX 

■L3*L3MAX 

•L4-L4MAX 


SIC  THE  PROGRAM  EXECUTES  OVER  THE  RANGE  OF  A DO  LOOP 
82C  FROM  E*0  TO  E*C7. 

83C  THE  VARIABLE  S SPECIFIES  THE  AVERAGE  SKIN  THICKNESS* 

84C  THE  VARIABLE  0 SPECIFIES  THE  RESISTIVITY  IN  0HM-CM  FOR  THE 
35C  TYPE  OF  MATERIAL  WHICH  COMPRISES  THE  SKIN* 

86C  FOR  EACH  ITERATION  A COMPUTATION  IS  MADE  OF  THE 
87C  FLUX  DENSITY  THE  TRANSFER  INDUCTANCE.  ANDTHE 
88C  TRANSFER  RESISTANCE. 

89 C ADDITIONALLY  .FOR  A SPECIFIED  LIGHTNING  WAVESHAPE 
90C  A TABULATION  OF  OPEN  CIRCUIT  VOLTAGE  VS*  TIME  IS  MADE* 

9IC  FOR  A TIME  PERIOD  T8  TO  T9  IN  STEPS  OF  T8  ( USECS) * 

92C  THE  WAVESHAPE  IS  CHARACTERIZED  BY  A DOUBLE  EXPONENTIAL  OR  A DAMPED 
93C  SINWAVE  MODIFIED  BY  THE  DIFFUSION  TIME  CONSTANT* 

94C  THE  CHOICE  OF  WAVESHAPE  IS  MADE  BY  THE  USER* FOR  A SINWAVE 
95C  SET  THE  VARIABLE  G4>I*0.ANY  OTHER  VALUE  DEFAULTS  TO  THE  DOUBLE 
96C  EXPONENTIAL *BOTH  TYPES  OF  EOUATIONS  ARE  SPECIFIED  WHEN  THE 
97C  USER  SELECTS  VALUES  FOR  THE  VARIABLES  "I  4".  "Gl".  "G2".  AND  "63". 

98C  SEE  THE  USERS  MANUAL  FOR  SUGGESTED  VALUES  TO  BE  USED* 

99C  FOR  WING--HORIZ--VERT  DATA  . (LINES  33-38)  . R IS  THE 
100C  LEADING  EDGE  RADIUS.C1  IS  THE  FWD  TO  AFT  LENGTH  OF  THE 
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101C  MAIN  BOX  (STRAIGHT  3ECTI0N)*T  18  THE  PMO  T»  APT  LENGTH 
108C  OP  THE  PLA»-S  (WINS),  TAPERED  TRAILING  EME  (HORIZ  STAB) 

109C  BR  RUDDER  CVERT  STAB) 

110  INTEGER  A*A1*  Aft*C7*  D4*D5*  D7*  08*  E*I*  J*K*LH«M*Ml*04*05*07«0G 

111  REAL  A8«B0*B1*B8»B3«BS»B6#B8*B9*C*C1* D*XMATD*D1*D8*D3*D6 
11B  REAL  E8*E6*  E7*E8*E9*P0*P1*P8«  P3*  P4»  F5*P6*  P7»P8*P9*  81*08*63 

113  REAL  H1*H8*H3*XNATI*I  1«I8«I3*14»IS*XJ*K1*K8*K3*K4*LI*L8*L3#L4*XL 

114  REAL  XN*XMATN*N7«H6*XNATN*0«01*03#86#Pl*P8*00*OI*88* 03*04 
US  REAL  R«R1«R8*RS*S* S7*T*Tl*T3*T4*T5#T7*TS*T9*  V.XNATV*  W 

116  REAL  X*X1«X6*XS*X7*X8*X9«Y*YI*Y6*  Y3*YS*Y8*  Y9*Z8»Z3 

117  90  P8RMATC  V) 

118  DIMENSION  PILES!  6) 

119  PI L INANE  PILE3/-6R0UCH0-*-HARP0-* -CHICO-* -8UNH0-*”ZEPP0M* -MAXWELL-/ 
180  D1NENS10N  XMATtK  1 6*  1 6)*XMATV<  1 6*  1 )*XMATI<  1 6*  1 )*XNATD(  1 6#  4) 

iei  DIMENSION  XMATN( 16*1 6)* XHATJ< 16) 
tee  REMIND  **6R0UCH0M 
tea  ENDPILE  -8R0UCH0" 

184  REMIND  "HARPO- 

185  ENDPILE  -HARPO" 

186  REMIND  "CHICO*’ 

187  ENDPILE  -CHICO- 

188  REMIND  -6UHN0- 
1S1  ENDPILE  -6UNH0- 

160  REMIND  -ZEPPO- 

161  ENDPILE  -ZEPPO* 

170  P1*3«141S987 

180  P8>6*88318S3 

190  E8-8. 71888 

800  Nl*l 6 

810  810  READ!  -MAXWELL-*  90)  A 

880  IPCA.E0«0)60T0S000 

830  80  TO  (840* 1840. i860. 1880*  50001* A 

840  840  READ(-MAXMELL-*90)A1*R1*R8*X1*Y1#C7»XS*  D3#  D4»  DS.YS*  06*  D7*  DO 
850  READ!  -MAX  WELL-* 90)  S*  XL* 0* L3*  03*  04*  OS* L 4* 06*87*  08*  TO*  T9*  1 4*81*88*  830 
861  PRINT  86e 

868  868  P0RNATC1H  • 17X«-«*  DIFFUSION- -COUPLING* -IN— FUSELAGE**-) 

870  PRINT  878 

871  PRINT  878 

878  878  P0RMATC1H/) 

880  PRINT  878 

886  IEUP*C7+1 

887  DO  1880  IED1M«1*IEUP 

888  IE*IEDUM-1 

890  REMIND  -GROUCHO- 

891  ENDPILE  -6R0UCH0- 

300  REMIND  -HARPO- 

301  ENDPILE  -HARPO- 

310  REMIND  -CHICO- 

311  ENDPILE  -CHICO- 
380  REMIND  -8UMM0- 
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381  ENDFILE  "6UMM0" 

330  S7-S 

350  IF< IE*E0*C7)  GOTO  183C 

340  X9«X5 

370  Y*«Y5 

380  LI  *L3 

390  L8«L4 

400  I FCD4.LE*  IE)  6070480 
410  68T0440 

480  480  X9«XS+<D3*(IE-D4)> 

430  IF<  D5*LT«IE)60T6S70 
440  440  1 F<  D7«LE* I E)  GOTO  440 
450  64T8480 

460  440  Y9*Y5*CD4*C  I E-D7)) 

470  I F(  08  «LT»IE)  60T0590 
480  480  IF(04.LE«IE)60T0500 
490  60T8S80 

500  500  LI*L3*(03*(IE*04)) 

510  1F<0S.LT*IE>60T0610 
580  580  IF(07«LE«1E)66T6540 
530  60T0440 

540  540  L8aL4+(06*(I  E*07>) 

550  IF<08«LE«IE)66TO64O 
540  64 TO  4 30 

570  570  X9aX5*< D3*D5) 

580  60T0440 

590  590  Y9«Y5*< 06*06) 

400  6078460 

610  610  LI  ■ L3*<03*05> 

480  6070580 

430  630  L2aL4+<06*08) 

440  640  HI*X|-R| 

650  K1«Y1“RI 
660  H8-RI 
670  K8-Y1-RI 
680  H3*R8 

490  K3«R8 

700  H4«X  I *R8 
710  K4*R8 

780  IF(A«EQ«AI) 0070740 
730  60704880 
740  740  CONTINUE 
741 C REM 

748  REWIND  M6R0UCH0M 

743  ENDFILE  "6R0UCHB- 

744  REWIND  "HARPO" 

745  ENDFILE  "HARPO" 

**50  D6770I a1 »N1 

760  WRITE  <"GROUCHO"* 90)1*3 
770  770  CONTINUE 
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WO  C»84<Y1-<RDR2D»89<X1-8*R8>*P19<RDR8> 

790  D-C/Nl 
800  X«0 

809  IUP»NJ»1 

810  088  50 1 DUN  » 1 # I UP 

811  I-IOUH-l 

880  Y»<Y1-R1>-I*I 

830  IF<Y<LE<Re>60T0860 

840  WRITE  <"HARP0"#90>X#Y 

850  8 50  CONTINUE 

860  860  IF<X<EQ<XDG0T0930 

870  X«X1 

879  IUP-ND1 

880  D0980I DUM«1« I UP 

881  I -I  DUN- 1 

890  Y*R2+D*1 

900  IF<Y<6E<Y1-RDG0T0930 

910  WRITE  <"HARP0"#90>X«Y 

980  980  CONTINUE 

930  9 30  Y>0 

940  D0980J*1#N1 

950  X8*Re«D«J 

960  IF<X8<GE<Xt-R8)G0T0990 

970  WRITE  C"HARPO**» 90>X8» Y 

980  980  CONTINUE 

990  990  IF<Y<EQ<YD60T01060 

1000  Y»Y1 

1009  IUP»N1*1 

1010  DOl  0501  DUM»  1 » I UP 

1011  I«IDUM-1 

1020  X-CX1-RD-D9I 
1030  IF<X<LE<RD60T01060 
1040  WRITE  <**HARP0**#90)X»Y 
1050  1050  CONTINUE 

1060  1060  KUP«N1«1 

1061  D01180KDUM-DKUP 
1068  K-KDUM-1 

1070  T1»K*D/R8 

1080  X8-H3-R8*C0S<T1> 

1090  Y2»K3*R8*SIN(T1 ) 

1100  IF(X8<6E<R8>Q0T01 130 
1110  WRITE  <MHARP0"*90>X8»Y8 
1120  1180  CONTINUE 
1130  1130  CONTINUE 
1131C  REM  RESET 

1139  KUP»N1*1 

1140  D01800KDUM«1*KUP 

1141  K>KDUM-1 
1150  Tt«K*D/R8 

1160  X8«H4*R8*SIN<TD 
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1170  Y8«K4-R8*C0S< T1 > 

1180  IF<V8<GE<R8>  G4T81 810 
1190  WRITE  < "MRRPO"» 90  > X 2#  Y 8 
1800  1800  CONTINUE 
1810  1810  CONTINUE 
181 1C  REM  RESET 
1819  LUP*N1*1 

1880  001880L0UM«1#LUP 

1881  L«LDUM-1 

1830  T1«L*D/R1 

1840  X8*H1»R1*C0S<  T1 ) 

1850  Y8“K1*R1*SIN<  T1 ) 

1860  ir<X8*I.E*Hl>00T0l890 
1870  WRITE  <MHARF0N»90>X8»Y8 
1880  1880  CONTINUE 
1890  1890  CONTINUE 
189 1C  REM  RESET 
1899  NUP*N1*1 

1300  001  360MDUM-UHUP 

1301  N-MDUN-1 

1310  T1«H*0/R1 

1380  X8*H8“Rl*SIN<  T1 > 

1330  Y8«K8«R1*C0S(T1> 

1340  IF<Y8<LE<K8J80T01370 
1350  WRITE  <"HARP0"#90>X8#Y8 
1360  1360  CONTINUE 
1370  1370  CONTINUE 
137 1C  REM  RESET 
1380  IF<X9<LE<R8JGOT01 480 
1390  IF<X9<LE<R1JG0T01430 
1400  IF<X9<GE<H4J80T01 540 
1410  IF<X9<6E<H1JG0T01 550 
1480  1 480  IF( Y9  <LE<R2J GOTO  1 450 
1430  1430  I F<Y9  <6E<K1 J GOTOI 500 
1440  GOTO  1480 

1450  1450  IF<X9<EO<Y9JGOTOI 690 
1460  X8*X9 

1470  Y8*Y9 

1480  X7»<R8”S8RT<R8*8”<<R8-Y9J »8J  J J 

1490  GOTO 1760 
1500  1500  X8*X9 

1510  X7«Rl-<SQRT<<RIJf8~<<RI>Yl*Y9J  tgJ  JJ 
1580  Y8«Y9 

1530  GOTO 1760 

1540  1540  IF<Y9<LE<R8JG0T0I 600 

1550  1550  IF<Y9<6E<KI > GOTOI 440 

1560  X8*X9 

1570  X7«X1 

1580  Y8*Y9 

1590  GOTOI 760 
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1600  1600  X8-X9 

1610  X7«X 1 ♦ < SORTC R2»  2- * Rg- Y9> * 2) > -R2 
1620  Y8*Y9 

1630  GOTO I 760 
1640  1640  X8"X9 

1650  X7»Xl*<S0RT(CRl)*2-CRl-Yl»Y9>t2>>-Rl 
1660  Y8-Y9 

1670  GOTO  1760 

1680  1680  IF<X9«GE*X1/2>GOT0173O 

1690  1690  X8"X9 

1700  X 7*0 

1710  Y8-Y9 

1720  60T0176O 

1730  1730  X8*X9 

1740  X7»X1 

1750  Y8aY9 

1760  1760  CONTINUE 

1761 C REM 

1770  ASSIGN  1778  TO  SW2900 

1772  GO  TO  2740 

1778  1778  CONTINUE 

1780  IF* I E* 6T*0> G0T01800 

1790  ASSIGN  1798  TO  SW3290 

1792  GO  TO  2910 

1798  1798  CONTINUE 

1800  1800  ASSIGN  1808  TO  SW4I80 

1802  GO  TO  3300 

1808  1808  CONTINUE 

1810  ASSIGN  1818  TO  SW4810 

1812  GO  TO  4190 

1818  1818  CONTINUE 

1820  1820  CONTINUE 

1830  1830  G0T021O 

1841  1840  PRINT  1842 

1842  18  42  FORMAT*  |H  *20X» "**DI FFUSI ON— COUPLING— IN—  WING****) 

1850  GOTO  1 890 

1861  1860  PRINT  1862 

1862  1862  FORMAT*  1H  * 10X*"**DI  FFUSION--COUPLING-- IN- -HORIZONTAL 
18 63*“ “STABILIZER**") 

1870  GOTO  1890 

1881  1880  PRINT  1882 

1882  1882  FORMAT*  IH  * 1 1X*"DI FFUSI0N--COUPLING--  IN— VERTICAL 
1883*““STABILIZER**"> 

1890  1890  PRINT  272 
1900  PRINT  272 

1910  READ*  "MAXWELL"* 90) A2»R» Cl »T#  S#XL#  C7#X  5»  D3»  D4#  D5«  YS*  D6#  D7»  DO 
1920  READ*  "MAX  WELL"* 90) 0*L3*  03*  04*0  5*  L 4*  06*07*  08*  T8»  T9*  I 4*  61*  G2»  63*  G4 

1930  REWIND  "GR0UCHO" 

1931  ENDFILE  "GROUCHO" 

1940  REWIND  "HARP0" 

Figure  29.  DIFFUSION  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  7 of  17) 


62 


DtrusitN 


02/04/75 


1941  ENDTILC  "HARPO" 

1950  REMIND  "CHICO" 

1951  ENDTILC  "CHICO" 

1940  REMIND  "8UNM0" 

1941  ENDTILC  "6UNN0" 

1970  S7-S 

1979  IEUP-C7M 

1980  DO  2720  IEDUN-1#IEUP 

1981  I£«IEDUN-1 

1990  IT(IE»E0*C7)60T0  2730 

2000  X9-XS 

2010  Y9-V5 

2020  L1-L3 

2030  L2«L4 

2040  IT< D4*LE*I E) 60T02040 
2050  00 TO 2080 

2040  2040  X9>XS*<D3*<1E-D4)> 

2070  IT( D5»LT«IE) 80T02210 
2060  2080  1T(D7*LE*1E)60T02100 
2090  00 TO 2 120 

2100  2100  Y9>YS*(D4*<IE-D7>) 

2110  1T(D6*LT*1E>OOT02230 
2120  2120  IT(04*LE*I  E)  60T021 40 
2130  60 TO  21 40 

2140  2140  Ll«L3*<03*<IE-04>> 

2150  IT(05*LT*IE)  60T026S0 
2140  2140  IT(07.LE.IE)«OT02180 
2170  60 TO 2280 

2180  2180  L2»L4*<04*<IE~07>> 

2190  I T<  08  *LE*  1 E)  60T02280 

2200  00T02270 

2210  2210  X9>X5*(D3«DS) 

2220  00T02080 

2230  22  0 Y9*Y5+(D6*D6> 

2240  68142120 

2250  2250  Ll«  L3*<03*05> 

2240  60T02140 

2270  2270  L2-L4*<04«08> 

2280  2280  IT<A*CO*A2>60T02300 

2290  00 TO 48 20 

2300  2300  M»ATAN(R/T> 

2310  REMIND  "QROUCHO" 

2311  ENDTILC  "OROUCHO" 

2320  C«(P14R)+2*Cl+2*( S0RT(R*2+T*2)> 

2330  D-C/Nl 
2340  002340N»l#Nl 
2350  MR1TE  ( "OROUCHO"# 90 ) N« S 
2340  2340  CONTINUE 

2370  REMIND  "HARM" 

2371  ENDTILC  "HARM" 
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2379  1 UP*N1 ♦ 1 

2380  0825701 DUH>1«I UP 

2381  1-IDUM-l 

2390  X2«(-T>»<I*D*C0S<W>> 

2400  Y2*l*D*SIN< W> 

2410  Y3--Y2 

2420  IF(X2»6T«O>80T0244O 
2430  G0T0253O 
2440  2440  Y2«R 
2450  Y3—Y2 

2460  1 F(X2*  GT*C1 ) G0T02480 
2470  G0T02S3O 

2480  2480  IFCX2-C1 •GT«R>80T0258O 
2490  Y2«S0RT((Rtf)>(X2-Cl>*2> 

2500  Y3—Y2 

2510  1F<X2«GT.C!*R)G0T02S9O 
2520  G0T0253O 

2530  2530  1F<Y2.EO*O>00T0256O 
2538  WRITE(**HARP0**#9O>  X2»Y2 
2540  WRITE  <**HARP0*%90>X2#Y3 
2550  Q0T0257O 

2560  2560  WRITE  (MHARF0M«9O)X2* Y2 
2570  2570  CONTINUE 
2S80  2560  WRITE  <"HARP0"»90>CI«>R*0 
2590  2590  IF(X9.6T.O)Q0T0264O 
2600  X8*X9 

2610  X7**<  T*< ABS<Y9/( TANC W) ) ) ) ) 

2620  Y8-Y9 

2630  G0T0267O 

2640  2640  X8-X9 

2650  X7«C1*S0RTCR»2-<Y9>»2> 

2660  Y8«Y9 

2670  2670  ASSIGN  2678  TO  SW2900 

2672  GO  TO  2740 

2678  2678  CONTINUE 

2680  I Ft  I E« GT.O) G0T02700 

2690  ASSIGN  2698  TO  SW  3290 

2692  GO  TO  2910 

2698  2698  CONTINUE 

2700  2700  ASSIGN  2708  TO  SW4180 

2702  GO  TO  3300 

2708  2708  CONTINUE 

2710  ASSIGN  2718  TO  SW4810 

2712  GO  TO  4190 

2718  2718  CONTINUE 

2720  2720  CONTINUE 

2730  2730  6OT021O 

2740  2 7 40  REWIND  "GROUCH0" 

2750  REWIND  "HARPO" 

2760  REWIND  "CHICO** 
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2761  EN0F1LE  "CHICO" 

2770  REWIND  "6UNM0" 

2771  ENDFILE  "6UMN0" 

2780  D92890I«1«N1 

2790  READ!  "6R9UCH9"* 90 >N»S 
2800  READ!  "HARP9"#90>X» Y 
2820  WRITE  <”CHIC9"»90>N«X»Y» S 
2830  Dl*SQRT<<X-X8)t2*<Y-Y8>t2> 

2840  D2*SQRT<  <X~X7)  *8*< Y*Y8>  *2) 

2850  IF<  D2* EQ«0) 69T92870 

2860  G9T928B0 

2870  2870  D2*S/2 

2880  2880  WRITE  <"6UMM0"*90>D1»  02 

2890  2890  C8NT1NUE 

2900  68  T8  SW2900 

2910  2910  REWIND  "CHIC9" 

2920  REWIND  "6UNM9" 

2930  REWIND  "6R8UCH8" 

2931  ENDFILE  "6R0UCH0" 

29  40  REWIND  "HARM** 

2941  ENDFILE  "HARP0" 

2950  REWIND  "ZEPP8" 

2951  ENDFILE  "ZEPPi" 

2964  READ<"CHIC9"«90>  < <XMATD<  IR9W*  I C9L)»  IC9L"1*  4)#IR8W»1#N1) 

2970  08  2978  IR8W>1«N1 

2974  WRITE! "GR8UCH9"»90>  tXMATD! IROW* 1C0L)» 1C9L" 1*  4) 

2978  2978  C9NTINUE 

2980  CALL  MATZER<XMATV#N1« 11 

2990  D93010I«1#NI 

3000  XHATV< I# l)«l 

»10  3010  CONTINUE 

3020  CALL  MATZER<XMATM»N1 »Nt) 

3030  CALL  MATZER<XMATI*N1*1> 

3040  CALL  MATZER<XMATN*N1»N1> 

3050  P2-6.28318 

3060  E2>2«  71828 

3070  RS» 30000 

307  6 307  6 F9RMAT<<4<1H  «G13*5>/>> 

3079  3079  F9RMAT<<5<1H  «Gl3*5>/>> 

3080  D93160I>1*N1 
3090  D931S0J>I»N! 

3100  IF<I *EC»J) 68T931 40 

31 1 0 R3»SQRT< < XMATD< I » C) -XMATD<  J#  2) > » 2* < XMATD< I » 3) -XMATD! J» 3>  > » 2) 
3120  XMATM<  I # J)  "ALOGt  R5/R3) 

3130  69T931  SO 

3140  3140  XMATM<I» J) ■AL96<RS*2/XMATD<  I#  4!)  > 

3150  3150  CONTINUE 
3160  3160  CONTINUE 
3170  CALL  MATI NV<  XMATM*  XMATN#  N 1 # N 1 ) 

3172  D9  3174  I9*l»Nl 
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31.3  XMATI C19# l>«0 

3174  3174  CONTINUE 

3175  DO  3179  I9>I#N| 

3176  00  3178  J9«I»NI 

3177  XMATI C 1 9 * 1 » »XMATI (19*1) *XMATNC 1 9 # J9 > 

3178  3178  CONTINUE 

3179  3179  CONTINUE 

3189  3189  F0RMATCSI3»5/) 

3190  11*0 

3eOO  DO 3220 I* 1#N| 

3210  1 1 »XMATI < I # I > + 1 1 

3220  3220  CONTINUE 
3230  D03280I>I#N| 

3240  READC  "GUMM0"#90>D1#  D2 

3250  XMATJC I >"XNATI C I# 1 )/I I 

3260  WRITE  C**HARP0"#90>  I#  Dl#  D2#XMATJC  I ) 

3270  WRITE  C"ZEPP0"#90)XMATJCI> 

3280  3280  CONTINUE 

3290  GO  TO  SW3290 

3300  3300  REWIND  "CHICO** 

3310  REWIND  "HARPS" 

3320  REWIND  "6UMM0" 

3330  REWIND  "ZEPPO" 

3340  D03830I>1*N1 

3350  READ!  "CHIC0"«90)XDUMN«X#Y# S 

3360  IF(I E* E0«0> GOTO 3400 

3370  READC  "GUMM0"#90)DI*D2 

3380  READC  "ZEPPO"# 90) XJ 

3390  G0T034I0 

3400  3400  READC  **HARP0"#93)II#D1#D2#XJ 
3410  3410  B8-CLl)/CDI4CSQRTCCLlt2)4CDl»2>>>> 

3420  B9«C  CXL-L1 ) )/C  D1*C  SORTC  CXL-Ll ) »2»C D! »2>  >) > 

3430  B1»CC 1 E* 5)4XJ)*C B8+B9) 

3435  A4«CL2”L1)*CABSCX7-X8)> 

3440  IFCX*EQ*X9)G0T03640 
3450  IFCY • E0*Y9) GOTO 3720 
3460  IFCX*LT*X9)GOT03560 
3470  IFCY*LT*Y9) 60T03S20 
3480  T3«ATANCCY-Y9)/CX-X9)> 

3490  Z2»-C  SINCT3)) 

3500  Z3»C0SCT3> 

3510  G0T0  3780 

3520  3520  T3-ATANC C Y9-YJ/CX-X9) > 

3530  Z2-SINCT3) 

3540  Z3-C0SCT3) 

3550  Q0T03780 

3560  3560  I FC  Y *LT*Y9)  G0T03610 
3570  T3«ATANCCY-Y9)/CX9-X)> 

3580  Z2*“C  5INC  T3> ) 

3590  Z3--CC0SCT3)) 
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•waa  OililTtA 

3410  3410  T3»ATANC<Y9-Y)/CX9-X>> 

3420  Z2*SIN( T3> 

3430  Z3*-CC0S(T3>> 

343S  00  TO  3700 
3440  3440  T3-PI/2 

3490  I F<  Y »LT* Y9> 00T03490 
3440  Z2*~(SIN<  T3>> 

3470  Z3-0 

3480  00T03780 

3490  3490  Z2*SIN<T3> 

3700  Z3*0 

3710  00T03780 

3720  3720  I F(X*OT*X9)OOT83740 

3730  Z2*0 

3740  Z3**l 

3750  02T03700 

3740  3740  Z2*0 

3770  Z3-I 

3780  3780  B2*B1«Z2 

3790  B3*B1*Z3 

3800  BS*B5+B2 

3810  84*B4*B3 

3820  BO*SORT<<BSt2)+<B4»2>> 

3830  3830  CONTINUE 

3840  I F( B5* E0 • 0) 88T83970 

3850  T 4*  ATANt ABS( 84) /ABS< BS) ) 

3840  IF<BS.0T.O)O8T83920 
3870  IF(B4.QT.0)68T83900 
3880  T5*180*<T4»S7«2958> 

3890  68 TO 40 10 

3900  3900  T5*180-<T4*57.2950> 

3910  08 TO  4010 

3920  3920  1F«B4.8T*0)68T039S0 
3930  T5*340"<T44  57*29  58) 

3940  60 TO 4010 

3950  3950  TS*T4»57#2956 

3940  68 TO 4010 

3970  3970  IF<B4*8T.0)68T84000 
3980  T 5*270 
3990  68 Tl  4010 
4000  4000  T5*90 

4011  4010  PRINT  4012 

4012  4012  F8RNATC  1H  »"HA6NETIC FIELD 

40134 CONFUTATION") 

4030  PRINT  272 

4031  PRINT  4032# X9 

4032  4032  FORMAT!  IN  #"X-C00RDINATE*"#613#4> 

4033  PRINT  4034#Y9 

4034  4034  FORMAT! IH  # "Y- COORDINATE*"# 61 3. 6) 
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4036  PRINT  40  37.L1 

4037  4037  FORMAT! 1H  *MZ1-C68RQINATE"M*G13'S) 

4036  PRINT  4039*1.8 

4039  4039  F6RNATC 1H  *MZ8~C66RD1NATE"M* 61 3. 5> 

4041  PRINT  4048 

4048  4048  FORMAT!  1H  *1X*"L88P  AREA  B-X  B-Y 

4043 A B-T8TAL  ANGLE**) 

4041  PRINT  4048 

4048  4068  FORMAT! 1H  # 39X# **<  WEBCRS/NE1  ER» 8)  (DEGREES)*4) 

4070  PRINT  878 

4081  PRINT  4088*A4*B5*B6*B0*T5 
4088  4088  FORMAT! 1H  * 5!  1H  *613*4)) 

4180  BO-O 

4130  B1>0 

41 40  B8«0 

41  SO  B3a0 

4160  B5>0 

41 70  B6a0 

4180  GO  T6  SW4180 

4190  4190  REWIND  **HARP6** 

4800  REWIND  **6UHH0N 
4810  REWIND  **ZEPP6** 

4880  XN«0 

4830  D84490I>1*N1 

4840  IF! IE«E0*0)G8T04880 

48  SO  READ(  **6UNM0**»90)  Dl*  D2 

4860  READ(  "ZEPP0M*90)XJ 

4870  G8T84890 

4880  4280  READ!  **HARP0***90>  II • 01*  D8*XJ 
4890  4890  04*! 1 E*9)*XJ 

4300  FI *S0RT!L8* 8+08*8) 

4310  F8sS0RT(L 1 *8+02*8) 

4380  F3-CL1-XL) 

4330  F4>!L8-XL) 

4340  F5*!F1-!L8>*!AL0G!!F1+L8>/D8>>> 

4350  F6"!F8-Ll*!AL08!!Fe+Ll)/D8>)) 

4360  F7«! SORT!  F3» 8+D8* 8) -F3*!  AL66C < F3+  SORT! F3»8+D8*8) )/D8) ) ) 

4370  F8"<  SORTC  F4»8+D8» 8)  -F4*(  AL86(  ( F 4+ SORT!  F4*8+ D8*8)  ) /D8>  ) ) 

4380  F9*SQRT!L2*2+D1 *8) 

4390  FO*SORT!Ll »8+Dl»8> 

4400  00*!F9-L8*AL0G!!F9+L8)/D1>) 

4410  O1>!FO-L1*AL0G!!FO+L1>/D1>) 

4480  08*! SORT!  F3»8+ Dl » 8)  -F3*!  AL66(  ( F 3+ SORT!  F3»8+ Dl  *8)  )/Dl ) ) > 

4430  03*!  SORTC  F4»8+D1 » 8) -F4M  AL6G(  < F4+S0RT!  F4»8+D1 » 8)  )/Dl ) ) ) 

4440  M7*F5-F6+F7-F8 

4450  N8-00-Q1+Q8-03 

4440  M7*M7*04 

4470  M8*M8*04 

4480  XM*!XM+!M7-M8)> 

4490  4490  CBNTINUE 
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4500  A8-S7*C 

4510  §1«0*<XL/A8) 

4520  PRINT  272 
4530  PRINT  272 

4541  PRINT  4542 

4542  45  42  FORMATdH  # 2X# "TRANSFER^ -♦-FUNCTION 
4543*^-* -♦-♦-♦-♦-♦-♦-♦- COMP UTATI ON") 

4560  PRINT  272 

4571  PRINT  4572 

4572  4572  FORMATdH  »8X* "TRANSFER  INDUCTANCE 

4573*  TRANSFER  RESISTANCE**) 

4581  PRINT  4582 

4582  4582  F0RMAT(1H  » 1 3X#  **<HENRI  ES)  C OHMS) **) 

4£9 1 PRINT  4S92«XM»01 

4592  4592  F0RNAT< 1H  » 12X# 61 3. 6» 22X# 613. 6) 

4601  PRINT  4602 

4602  4602  F0RMATMH  #**OPEN  CIRCUIT  VOLTAGE**) 

4610  PRINT  272 

4ol  3 PRINT  4614 

4614  461  4 FORMATdH  «**TIME  VOLTS**) 

4619  T7-0 

4620  DO  4720  1 DUMMY-1 *999 

4621  T7-T7> T8 

4622  IF  <T7.GT.T9)  GO  TO  4721 

4623  IF  (G4.E0*1.0)  GO  TO  4692 

4630  12*1 4*<  <-Gl*EXP<-Gl*T7)  ) ♦ 62*EXP(  -62*T7)  ) 

4640  12-12*1 4* < <G1*G3)*EXP<<-G1“G3)4T7)) 

4650  1 3-1 4*<  G2*G3)*£XP<  <-G2“G3)*T7) 

4660  1 5-12*13 

4670  E7-01  *1  4*<  EXP< -Gl  *T7)  - EXPC  - G2*T7)  )*<  l-EXPC  -G3*T7)  ) 

4680  E8-XM41 5 

4690  EV-E7-E8 

4691  GO  TO  4711 

4692  4692  AMP-(  I 4*SIN<2.-P  1 -G1-T7)  )*EXP<  -G2-T7)  *<  1-EXP(  - G3*T7)  ) 

4693  DAMPDT-I  4*<2.*P1*G1*C05(2«*PI*GI*T7>  ) 

4694  DAMP DT- DAMP DT*< EXPC -G2*T7>- EXPC <-G2-G3)4T7)  ) 

4695  DAMP-I4*<SIN<2.*P|*GI*T7)) 

4696  DAMPT-DAMP4CC  62* G3)* EXPC  < -G2-63)*T7)  -61*EXP<  -GI*T7)  ) 

4697  DAMP DT- DAMP DT* DAMP T 

4698  E7-014AMP 

4699  E8-XM*DAMPDT 

4700  E9-E7-E8 

4711  4711  PRINT  4712#T7»E9 

4712  4712  FORMATdH  .G13.6.3H  .613.6) 

4720  4720  CONTINUE 

4721  4721  CONTINUE 
4730  PRINT  272 
4741  PRINT  4742 

47  42  4742  FORMATdH  .75C1H-)) 

4760  PRINT  272 
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PRINT  878 
PRINT  878 
XH-0 
• l»0 

60  T0  SW40IO 
4080  PRINT  878 
PRINT  40  38 

4038  FORMAT! 1H  » "DATA  READ  STATEMENT  O0ES  N0T  CONTAIN") 
PRINT  4048 

4048  P0RNATC 1H  ."VALUES  WHICH  CORRESPOND  TO  THIS**) 

PRINT  40S8 

4058  FORMAT! 1H  . **GEOHETRY • CHECK  ALL  DATA  STATEMENTS**) 
PRINT  4068 

4068  FORHATC 1H  . "TO  BE  SURE  THAT  THEY  ARE  CONSISTENT**) 
PRINT  4078 

4078  FORMAT! 1H  ."WITH  HIE  GEOMETRY  YOU  ARE  EVALUATING.**) 
5000  STOP) END 

SUBROUTINE  MATZERIXM. I ROW# ICOL) 

DIMENSION  XMIIROW.IC0L) 

DO  3800  19* 1. IROW 

DO  3190  J9>1*IC0L 

XMII9.J9)*0 

3190  CONTINUE 

3800  CONTINUE 

RETURN 

END 

SUBROUTINE  MATINVI  A.  B*  IROW*  1 COL) 

DIMENSION  AIIROW.  ICOD.BI  IROW*  ICOL). Cl  1*1) 

DO  3540  19*1  .IROW 
DO  3530  J9«l  . ICOL 
BII9.J9)«AII9.J9> 

3530  CONTINUE 
3540  CONTINUE 

CALL  MATRIX! 6. B. A*  C. 1 RO W. I COL. I RO W. I COL. ICOL) 

RETURN 

END 

SUBROUTI NE  MATRI X I I OP. A. B»  C. 1 » J. K. L. M) 

REAL  A.B.C.TEMP 

DIMENSION  AII.J).B!1*J).CII«J) 

DIMENSION  LABEL! 16) 

GO  TO  ! 101. 108. 103.104. 800. 300. 4C0).  IOP 


101  ASSIGN  111  TO  IP) 
108  ASSIGN  118  TO  IP) 


GO  TO 
GO  TO 
GO  TO 


103  ASSIGN  113  TO  IP)  GO  TO  100 

104  ASSIGN  114  TO  IP 
100  DO  180  1 1* l.K 

DO  180  18*1. L 

GO  TO  IP.!U1. 118.113*114) 

111  C! I 1* I8)"A! Il»18)+B! 11*18))  GO  TO 
118  C!I 1 * I8)»A! II. I 8) -Bill. 18)1  60  TO 
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6140  113  C<I I# 1£)*A<I 1# 1£>*8<1 1# 1£>)  GO  TO  l£0 
6150  114  C<  1 1* 1£)*A< I 1*  1 £>/£(!  1*  1£> 

6160  l£0  CONTINUE 
6170  GO  TO  500 
6180  £00  DO  £10  ll»l*K 
6190  00  £10  I£*1»L 
<£00  TEMP»0< 

<£10  DO  £05  13*1»M 

<££0  £05  TEMP»TEMP+A<  1 1#  1 3)*B<  13#  12) 

<£30  £10  C<11#1£)»TEMP 
<£40  GO  TO  500 
<£50  300  NR»Kf  NC=K 
<£60  DO  21  J 1 » 1 # NR 
<£70  £1  LABEL<Jl)>JI 
<£80  00  291  J 1 • 1 # NR 
<£90  TMP1-0. 

6300  00  121  J2*J 1#NR 
6310  TNP2*CABS( A<  J£#  J 1 ) ) 

6320  IFCTHP2-TNP1)  121# 121* 1210 

6330  1210  THP1-TMP2 

6340  IB1  G»J2 

6350  1 £ 1 CONTINUE 

6360  IF<IBIG<EO<JD  60  TO  201 

6370  00  141  J2-DNC 

6380  TEMP>A<Jl#J2) 

6390  A<J1«J2)>A<1B1G#J£> 

6400  1 41  A<IBI  6#  J2>»TEHP 

6410  I>LABEL<JD 

6420  LABEL<JD*LABEL<1BI G) 

6430  LABEL<IBIG)*I 

6440  201  TEMP«A<J1#JD 

6450  A<J1#JD»1<0 

6460  DO  2£1  J2-DNC 

6470  221  A< Jl# J£)<A< Jl# J21/TEMP 

6480  DO  281  J2*DNR 

6490  1F<J2<E0<J 1)  GO  TO  281 

6500  TEMP*A<J£*JD 

6510  A<J£#JD«0< 

6520  DO  241  J3<  DNC 

6530  £41  A<J£#J3XA<J£#J3)-TEMP*A<J1»J3> 

6540  281  CONTINUE 

6550  £91  CONTINUE 

6560  301  N1 *NR” 1 

6570  DO  391  J1«1#N1 

6580  DO  3£1  J£«JDNR 

65S0  1F<LABEL<J21<NE<J1 > GO  TO  3£1 

6600  IF<J2<E0<JD  GO  TO  391 

6610  GO  TO  341 

6620  321  CONTINUE 

6630  341  DO  361  J3*1»NR 
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6640  TEMP*A(J3»J1> 

6650  ACJ3»J])*A(J3,J2> 

6660  361  A(J3»J2)«TEMP 
6670  LA8EL(J2)*LA8EL( Jl) 
6680  391  CONTINUE 
6690  6f  T8  500 
6700  400  M 410  1 1*1#K 
6710  M 410  12>1«L 
6720  410  C(12»l 1>*ACI 1*12) 
6730  500  RETURN! END 
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etries.  Once  the  particular  geometry  being  evaluated  has  been  defined,  the 
program  defines  a horizontal  plane  containing  both  the  electrical  circuit  con- 
ductor and  a structural  return  path,  and  stores  the  coordinates  of  this  plane 
for  utilization  later  in  the  program.  This  is  done  in  lines  1380-1760  or  2590- 
2708,  At  this  point  in  the  program,  the  information  that  has  been  generated 
is: 


• Location  of  the  enclosed  electrical  circuit  conductor 

• Coordinates  of  all  current  carrying  filaments  representing  the  geo- 
metrical structure 

• Loop  or  surface  for  which  flux  linkage  is  to  be  computed 

The  program  now  branches  to  a subroutine  (lines  2680-2900)  which  determines 
the  distance,  Dl,  from  each  of  the  skin  current  filaments  to  the  enclosed  con- 
ductor, and  additionally  computes  the  distance,  D2,  from  each  of  the  skin  cur- 
rent filaments  to  the  return  skin  conductor.  These  values  are  retained  and 
stored  in  one  of  the  files  for  later  use. 

The  next  computation  to  be  performed  is  the  current  distribution  at  each 
filament.  The  computer  program  partitions  the  input  lightning  current  and 
defines  for  each  of  the  current  filaments  a fractional  portion  of  the  total  cur- 
rent. The  portion  of  the  lightning  current  assigned  to  each  current  filament 
depends  on  the  geometry  and  the  location  of  the  current  filament  in  that  geom- 
etry; this  operation  is  executed  in  lines  2910-3290.  At  this  point  in  the  pro- 
gram execution,  the  computer  program  has  defined  the  location  of  each  current 
carrying  filament,  the  current  distribution  in  each  of  these  current  filaments, 
and  the  location  of  the  enclosed  electrical  circuit  conductor.  It  has  also  defined 
a loop  through  which  flux  linkages  are  to  be  computed. 

The  program  now  branches  to  the  subroutine  in  which  flux  density  equa- 
tions an  used  to  compute  the  flux  density,  its  vector  components,  and  its 
orientation  at  a given  point  inside  the  geometry  of  interest.  The  point  se- 
lected for  this  computation  depends  upon  the  user's  selection  of  the  enclosed 
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circuit  conductor  initial  location,  Zl.  The  computational  operations  to  obtain 
the  flux  density  are  performed  in  lines  3300-4 1'80. 

The  computer  program  then  branches  to  the  subroutine  that  computes  the 
transfer  inductance,  M,  and  the  transfer  resistance,  R.  In  execution  of  this 
subroutine  a computation  is  made  of  the  open  circuit  voltage  versus  time,  util- 
izing the  transfer  functions.  Computation  of  the  transfer  inductance  is  per- 
formed by  reading  in  the  previously  computed  and  stored  values  of  D1  and  D2 
as  well  as  the  value  of  the  current  distribution  for  each  of  the  current  filaments. 
These  values  are  inserted  into  the  flux  equation  for  each  filament  out  to  a dis- 
tance Dl,  which  is  subtracted  from  the  flux  that  is  computed  for  the  skin  cur- 
rent filament  out  to  the  distance  D2.  This  difference  in  flux  is  the  net  flux 
linking  the  defined  plane.  The  transfer  inductance  is  the  summation  of  all  of 
these  individual  fluxes  divided  by  the  total  lightning  current  that  flows  through 
the  complete  structure.  These  operations  are  performed  in  lines  4190-4490. 

The  transfer  resistance  is  computed  using  the  equation 

R = pL/A  (63) 

where: 

p = resistivity  of  the  skin  material  (ohm-cm) 

A = cross-sectional  area  of  the  geometry  skin  (ma) 

L = total  length  of  the  geometry  being  evaluated  (meters) 

This  is  obtained  from  lines  4500-4592.  After  computing  the  transfer  induc- 
tance, M,  and  the  transfer  resistance,  R,  the  program  computes  the  induced 
voltage  in  the  specified  electrical  circuit,  utilizing  a lightning  waveform 
described  by  the  user's  dat^  inputs.  Since  naturally  occurring  lightning  strokes 
vary  greatly  in  waveform,  the  user  may  select  the  waveform  for  which  pro- 
tection is  to  be  designed  (one  for  either  a damped  oscillatory  or  a double  ex- 
ponential waveform  may  be  used).  The  waveform  of  the  portion  of  lightning 
current  appearing  at  the  inside  surface  of  the  skin,  and  thus  in  the  skin  cur- 
rent filaments  described  by  this  program,  is  not  the  original  lightning  current 
waveform.  Instead,  it  is  modified  by  a diffusion  time  constant  in  the  manner 
described  in  Reference  5.  This  is  accomplished  directly  in  the  induced  vol- 
tage equation.  The  resulting  open  circuit  voltages  are  then  computed  and  tab- 
ulated as  a function  of  time.  The  user  has  control  over  the  time  duration  and 
increments  over  which  this  computation  is  executed.  These  operations  are 
performed  in  lines  4601-4810. 

After  completing  this  computation,  the  program  loops  bark  and  determines 
if  modifications  to  the  previous  data  set  have  been  requested.  If  modifications 
are  to  be  made,  program  execution  repeats,  using  this  new  data  set.  A new 
set  of  transfer  functions  is  then  computed  along  with  the  corresponding  open 

circuit  voltages. 

If  no  modifications  to  the  data  have  been  requested,  the  program  then 
determines  if  another  geometry  has  been  selected  for  evaluation.  If  such  a 
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geometry  is  to  be  evaluated,  the  constants  of  that  geometry,  the  initial  location 
of  & circuit  conductor  inside  that  geometry,  and  the  modifiers  which  will  be 
used  to  reposition  or  relocate  that  electrical  circuit  conductor  are  read  and 
the  program  operates  as  before.  After  all  modifications  in  all  geometrical 
models  have  been  completed,  the  program  reaches  an  end. 

The  output  data  returned  from  the  program  are  the  coordinates  of  the 
circuit  conductor  for  the  execution  in  progress,  flux  density,  flux  density 
orientation  and  vector  components,  transfer  inductance,  transfer  resistance, 
and  a table  of  lightning  induced  voltage  versus  time  in  the  open  circuits  of 
interest. 

VALDATION  OF  DIFFUSION 


The  criteria  used  to  evaluate  the  validity  of  the  computer  program  were 
(1)  to  determine  if  it  returned  the  same  answer  that  could  be  manually  com- 
puted for  a textbook  calculable  geometry,  and  (2)  to  compare  computer  gen- 
erated values  to  those  of  aircraft  on  which  experimental  measurements  are 
available  from  which  empirically  derived  transfer  functions  were  available. 
Two  illustrative  cases  were  selected  and  are  presented. 

CASE  1 

The  object  in  to  compute  the  mutual  inductance  between  a single  current 
filament  and  a loop  with  a configuration  as  shown  in  Figure  30. 


Figure  30.  Single  Current  Carrying  Filament  and  Circuit  Loop 

The  expression  which  determines  the  flux  linking  the  loop  is  obtained 
for  this  case  from 


or 


(64) 


♦ is 


in 


b + D 
D 


(65) 
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The  mutual  inductance  is  obtained  by  dividing  by  1;  thus, 

t _ ^la  _ ^ oa  •-  b + D 
Li a " - « - 


2*  ^ a 


(66) 


Values  were  selected  for  this  geometry  as  follows: 


a = 50  cm;  b = 217  cm;  D = 150  cm 


and  Lia  was  computed  to  be  8.  9 x 10“®  henries. 


The  computer  generated  value  for  this  case  (7.2  x 10“®  henries)  is  presented 
in  Figure  31. 

CASE  2 


The  object  is  to  compute  the  flux  linking  an  electrical  circuit  that  is 
centered  in  the  cylindrical  fuselage  as  shown  in  Figure  32.  Because  of  the 
symmetry,  the  total  flux  linking  this  plane  should  be  equal  to  zero.  The 
computer  generated  results  are  shown  in  Figure  33. 


It  is  evident  that  in  the  limit,  as  the  modeled  geometry  is  simplified,  the 
analytical  expressions  evaluated  by  the  computer  program  DIFFUSION  reduce 
to  easily  computed,  classical  formulas. 
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Figure  31.  Diffusion  Calculated  Values  of  Transfer  Functions 

and  Open  Circuit  Induced  Voltage  in  Single  Geometry 
of  Figure  30 


Figure  32.  Electrical  Circuit  Loop  Inside  a Cylindrical  Fuselage 
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Section  3 

APERTURE  FIELDS 


APERTURE  THEORY 

EQUIVALENT  MAGNETIC  DIPOLES 


If  a magnetic  field  exists  tangentially  to  a surface  in  which  an  aperture 
exists,  the  fields  induced  on  the  other  side  of  that  aperture  may  be  treated 
those  induced  by  a dipole  of  appropriate  strength  lying  in  the  plane  of  that 
aperture  (Figure  34).  A mathematically  tractable  aperture  is  an  ellipse  of 


External  Area 


External 

Magnetic 

Field 


Y 

4 


♦ z 


j ► Internal  Area 

I 

I 


Equivalent  Dipole  Applicable 
to  Internal  Area 


Figure  34.  Equivalent  Dipole  Presented  by  Aperture 
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Preceding  pege  blank 
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The  alv  and  are 
the  elliptical  aperture  shown 


factors  related  to  the  shape  of  the  aperture, 
in  Figure  35: 


2ir 

(i1/8)sea 

(68) 

3 

K(e8)  - E(ea) 

2w 

(f  1/2)®ea(l  - ea) 

(69) 

3 

E(eV-  ( 1 - ea)  K(ea) 

2 Tt 

(ij/^O-e8) 

(70) 

3 

E(ea) 

For 


where: 


and  K(ea)  and  E(ea)  are  elliptic  integrals  of  the  first  and  second  kinds, 
respectively. 


At  the  surface  of  a conductor  the  Z component  of  H(ext),  Hz(ext)  must 
vanish  if  either  the  conductance  is  high  enough  or  ti,e  frequency  of  concern 
is  high  enough  so  the  skin  depth  is  small  compared  to  the  thickness  of  the 
conductor.  For  the  cases  of  present  interest  component  Hzwill  frequently  be 
zero,  by  virtue  of  the  geometry  of  the  current  flow  producing  magnetic  field 
H(ext).  Likewise,  for  the  cases  of  present  interest,  the  frequencies  at  which 
a magnetic  field  may  penetrate  in  a Z direction  are  low  enough  that  they  are 
not  of  concern.  Accordingly,  assume  that  Hz(ext)  = 0.  Under  these  conditions 
the  equivalent  dipoles  are: 


• Equivalent  dipole  lying  along  the  X axis: 

M,  = (H£),  = onHx(ext)  (71) 

• Equivalent  dipole  lying  along  the  Y axis: 

My  = (Hf)y  - ossHy(ext)  (72) 


For  the  case  in  which  the  major  axis  of  the  elliptical  aperture  is  oriented 
along  the  X axis  (as  shown  in  Figure  35),  the  values  of  an,  ox,  and  a 33  are  given 
on  Figure  36.  If  the  major  axis  is  oriented  along  the  Y axis,  the  same  curve 
is  applicable  if  the  designations  of  on  and  as s are  reversed. 


FIRST  ORDER  DIPOLE  APPROXIMATION  TO  INTERNAL 
MAGNETIC  FIELD 


Considering  a magnetic  dipole  of  strength  IIvi  located  along  the  Y axis, 
the  coordinate  geometry  would  be  as  shown  in  Figure  37.  At  point  P the  mag 
netic  potential,  M,  is. 
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(73) 


81 


0 


0.2  0.4  0.  r, 

1,/tJ 


0.8  1.0 


‘ Figure  36,  Shape  Factor  for  Elliptical  Apertures 
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M Details  of  Dipole  Located  at  Origin 


Figure  37.  Magnetic  Dipole  Oriented  Along  Y Axis 
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where: 


(For  purposes  of  clarity,  constant  Kj  is  reserved  for  a later  formulation  with 
the  dipole  located  along  the  X axis. ) 


Dipole  theory  generally  assumes  that  point  P is  sufficiently  far  from  the 
origin  that  ra  and  ra  may  be  approximated  (Figure  38)  as: 

t. 


ra  = r - -2-  cos0 


r9  = r + -g-  cos  0 
* 2 


Under  these  circumstances: 

M = K0 


r - -y-  COS0 


-H— ' 1 

r + ~2~  c°s0j 


M = K, 


If  S.J2  <Kr,  then: 


Jt  S. 

r + -£■  cos0  - r + -*■  cos 0 

Cd 

ra  -f  ^-!cos0  3 
^ I 


K,i,cos0 

M = ■ ■ 8 -s= 

r3 


(74) 

(75) 

(76) 


(77) 


(78) 


or: 


M 


(Hi,)cos0 
■ I;  r8" 


(79) 


Vx2+  z2 


Figure  38.  Approximations  Used  in  Elementary  Dipole  Analysis 
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HIGHER  ORDER  APPROXIMATION  TO  INTERNAL  MAGNETIC  FIELD 


i 

\ 

l 

§ 

l 

‘ 
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The  above  formulation  is  valid  when  the  point  at  which  the  field  is  to  be 
calculated  is  at  a distance,  from  the  aperture,  that  is  large  compared  to  the 
dimensions  of  the  aperture.  If  an  attempt  is  made  to  calculate  the  fields  close 
to  the  aperture,  the  results  are  inaccurate.  If  the  distance  from  the  point  to 
the  aperture  goes  to  zero,  the  fields  become  infinite,  whereas  they  can  never 
in  fact  become  larger  (barring  reflections)  than  the  external  field.  It  can  be 
shown,  in  fact,  that  the  actual  field  strength  in  the  plane  of  the  aperture  (as 
distinct  from  the  strength  of  the  equivalent  dipole)  will  be  half  the  field  strength 
that  would  exist  if  the  aperture  were  not  there. 

To  maintain  a little  more  numerical  accuracy  near  the  aperture,  the  ap- 
proximations made  by  Equations  74  and  75  will  not  be  made,  but  Equation  73  will 
instead  be  expanded  by  a power  series. 

DIPOLE  ORIENTED  ALONG  Y AXIS 


In  Figure  37: 


or: 

where: 

Ca 

Thus: 


X3  + z2  + |y  - 3 

(80) 

fx3  + Z3  + Jy  4-  ^-J  3 

(81) 

= <C8  - Yf,)1^ 

(82) 

= <CS  + Yfg)1'3 

(83) 

Xa  + Y3  + Za  + 


M = K3 


<CB  - Y l2)-''s  - (Ca  + Y £an'* 


(84) 


Expanding  by  the  binomial  theorem  and  combining  like  terms: 


M = K 


b„C“1/a  t-  b1(Yje?)C„-3/a  t-  b^Yf,)3^'3 

L + b3(Y£s)3C^'2  + 

" hoCa-1'3  + VYL,)Ca-*'a  - b?(Yis)aCiT6'a 
+ ^(Yi^C^3  +■ 


where: 
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h0  = 1 


L 


K4 


j 

j 


j 


I 


i 

j 


_1  3 _1_  3 

Da  ' 2*2*2!  "8 

_ 1 1 j>  _1_  _ _15  _ J5_ 

b3  = 2 * 2 ' 2 * 31  = 48  ” 16 

1 3 5 7 1 105 

b*  = 2 * 2 ’ 2 ’ 2 ‘ 4!  "384 

2 3 5 7 9 J_  _ 840  _ _63_ 

bs  = 2 ' 2 ' 2 ’ 2 ' 2 ' 5!  ' 3840  256 

1 3 5 7 9 U _L  _ 10395 

b«  = 2 * 2 ’ 2 ’ 2 ’ 2 ’ 2 ’6!  = 46080 

135135  3003 

b7  = 645120  14336 


Thus: 


: 

'•'9 

4 


(86) 


2M  = C«'a (b^a - Ylfl)  (3/2)  C ^ (2 Y - la) 


C 7'a(ba)  (2Yala)  - ba(Yia)a(7/2)Cae/a(2Y  - la) 


Ci4/8 


+ CU'a  (bR)  (5Y4ia5)  - bB(Yia)5(ll/2)  Ca9'a(2Y  - la) 


paa/a 

a 


+ Ca^a  (b7)  (7 Yeig7)  - b7(Yife)7  (15/2)  C^/a  (2Y  - l9) 


pao/a 

'“'3 


2Ka  bt 


3(Yla)(Y-la/2)  - Cya  ^ 7(Yia)»(Y-ia/2)  - 2CaYal, 


+ b - 5CaY4laS 


p la/a 


15(Yia)7  (Y  - ia/2)  - 7CaY6ia 


ca17/a 


The  gradient  in  tlie  Z direction  is: 


The  partial  differentiation  follows  the  same  format  as  Equation  22; 
therefore: 


„ , , . ^ 6 bjtYiglZ  ^ 14  ba(Yia)*Z 

rl " 8 1 — cp~  — cj7i — 


+ 22  bs(Yla)5Z  + 30b7(Yla)7Z  + 


cp >* 


Ca17'a 


DIPOLE  ORIENTED  ALONG  X AXIS 


Following  the  identical  line  of  attack,  with  the  dipole  oriented  along  the 
X axis,  yields  the  following  relationships: 

H | I , 2K,  rb  3<Xi,)(X-y2)-C,l,  + 7<X2,m-i,/2)-2C,X»i« 

Tl  1 [ 1 s C p 

+ b H(X^>MX-f,/2)-5C,X^  . . 

6 Cia/a 

+ h 15(XJ(!1)7(X-i1/2)  - 7C1X«i  7 + 

7 C117/a 
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, fe^tXI^Y  + 14  bjCXij^Y  + 22bB(X*1)BY 
Ll  [ c?5  Cj75  CpTa 


+ SOhfOU^Y  + 


C 17/a 


BbjOM^Z  + 14bsC 

cj73  Cjs 

+ SOb^OCjg^Z  + _ 


(Xl^Z  + 22b5(Xi1)EZ 


Ci3/a 


p i7/a 

^1 


TOTAL  MAGNETIC  FIELD 

The  total  field  at  point  P is  that  due  to  the  sum  of  the  external  fields  in 
the  Y and  X directions: 


Hx  = H*j,|+Hx,V 

Hv  = H./j+fM, 
If  \ > 
Hz  - HjL  I + Hj L] 


(100) 


[3x(x-^J  - CJ  + C^*K#«b[3YX] 

+ |ci^/aK1jf13  RX3|x-^-J  - 2C1Xa]+|ca^/aKai!3[7Yi*X] 


+"^crt3/aVi5[llxB  x_^ 


5C  jX4]  + C ?#»  £ [1 1 Y^X  J 


+ fy|fcr17/a!fV17  t.l5X7  [x-|i  - 7C1X6]  + |Yf|cip7/aK5rC7[l5Y7X] 


Cfe/9K1£l  Coxy]  + ca-*/aKaia  [sy  j y - ^ - ca] 

+ |c^/aK1£3[7X3Y]  + |cs-9/*Ka£a3[7Y3  (\  - 4r  I “ 2C aY33 


+ ^c-3/aK1iB[llxBY]+^cn3/aV3SBiYs(Y.^i  . 5CaY4] 

+ mlCi~17/3Ki£i7  [ISX’Yj  + l^lc^K^7  [15Y®  j Y j - 7CaYs] 
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(103) 


Hz  = C1~6/aK1J?1  [3XZ]  + [3YZ] 

5 

8'1  "i~i  u J 8 


+ [7X3Z j + ICa^K^a3  [7Y3Z] 


+ -^Ci'13/8Kiii6  Cl  1X6Z ] + C^13/a Kgig  [1 1Y6Z] 


128 

3003 


3003 


7 [l 5X7Z ] + ^Ca^K^a7  [15Y7Z] 

Equations  101  through  103  may  be  placed  in  a format  more  suitable  for 
machine  calculation,  as  follows: 


Hx  = GjXFj  + GjXFa-GjX  F3 
Hy  = G4  x Fj  + Gg  x Fa  - Ge  x F4 
Hz  = G7  x Fj  + Ga  x Fa 

where: 

G = 3Kt  (X  - i/2) 


P 1*5 


3KjX 


P 1.6 

'-a 


G.  = -iSlfi 

3 p 1.6 

g4  = *Eli 

Ci-e 


n = 3Ka  (Y  - ij2) 
6 Cg1*8 

Ga  = -Ml 

8 C 1.6 


g7=^£ 

7 C,1-8 


3KaZ 

C,1-6 


'.■feb- 


Yi  \ 

Fa  = —J  + 

« r- 1 


1.45833 
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Xij 

I cTI  1 

Ci  i 

. 45833  1 

Yj U 3 

+ 1.804688 

Yi,  V 

ca 

C ./ 

(104) 
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+ 2.  094727  _i 
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; Yi- V 9 7 
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Fa  = 1 + 2(Xll)a  + 5(Xf  / + 7(Xi1)«  . . . 

F4  = 1 + 2<Yi„)3  + 5(Yfa)4  + 7(Yj?8)8.  . . 

In  Figure  37,  product  Hi  is  the  strength  of  the  equivalent  dipole  produced 
in  the  aperture  by  the  external  magnetic  field.  For  the  portion  of  the  field 
caused  by  the  component  of  the  field  lying  along  the  Y axis,  the  strength  of 
the  dipole  is: 


(Hl)r  = o-aaHylext) 

(107) 

and  for  the  dipole  lying  along  the  X axis,  the  strength  is: 

(Hf)x  = o-jj  H„(ext) 

(108) 

The  dipole  moment  is  appropriate  for  use  in  the  classical  dipole  formula- 
tions based  on  Equation  86  but  is  not  appropriate  for  the  power  series  formu- 
lation used  in  the  main  text  of  ihis  report.  If  the  £ factor  is  taken  in  the  der- 
ivations as  the  actual  physical  dimension  of  the  aperture,  then  for  and  Ka: 

l'  _ hj c . °n  H„(ext) 
1 4 :r  j 

(109) 

and: 

Hy  anHv(ext) 
= = 

4ir  4w£a 

(110) 

where  fj  and  la  are  the  major  and  minor  dimensions,  respectively,  of  the 
elliptical  aperture. 


REFLECTING  SURFACES 
Two  Parallel  Plates 

The  problem  of  field  penetration  into  the  region  between  two  parallel  plates 
is  of  considerable  interest,  because  it  applies  to  the  degradation  of  shield  in- 
tegrity caused  by  the  presence  of  small  apertures.  The  preceding  analysis 
may  be  extended  to  two  parallel  plates,  one  having  an  aperture  and  the  other 
continuous,  by  using  image  theory. 

The  image  of  the  electric  dipole  moment  is  colinear  with  the  dipole  vector: 
however,  the  image  of  the  magnetic  dipole  is  antiparallel  with  the  magnetic 
dipole  vector.  Taking  this  into  consideration,  a doubly  infinite  array  of  images 
is  formed,  as  shown  in  Figure  39.  The  field  components  at  a particular  point 
in  space  may  be  obtained  by  an  algebraic  addition  of  all  of  the  contributions 
from  the  aperture  dipoles  and  the  image  dipoles. 

Multiple  Reflecting  Surfaces 

In  principle,  a rectangular  area  or  volume  behind  the  aperture  could  be 
formed  by  two  or  four  additional  reflecting  surfaces,  as  shown  in  Figure  40. 
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Figure  40.  Multiple  Reflecting  Surfaces 
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Although  the  number  of  images  increases  greatly,  only  the  images  fairly  near 
the  surface  generally  need  to  be  considered.  Multiple  reflecting  surfaces  have 
not  been  incorporated  into  this  program. 

FLUX  LINKING  A LOOP 

Figure  41  shows  a loop  defined  by  four  points  in  spaces  PI  through  P4,  all 
of  the  points  being  assumed  to  lie  in  the  indicated  plane.  At  some  point  PL 
within  this  loop  there  will  be  a magnetic  flux  vector,  T?.  The  XYZ  components 
of  this  vector  may  be  determined  from  the  previous  equations.  The  component 
of  H that  is  normal  to  the  plane  is  that  part  parallel  to  the  normal  vector,  TT, 
at  point  PL,  or: 

Hn  = H coso  (111) 


Figure  41.  Flux  Linking  Arbitrary  Four  Sided  Loop 


»i 


wnich  is  equivalent  in  vector  notation  to  the  dot  product: 

H„  = H # NlT  (112) 

where  NU  is  the  unit  vector  normal  to  the  plane  defined  by  points  PI  through  P4. 


Points  PI,  P2.  and  P3  will  be  used  to  define  the  plane  in  which  all  of  the 
points  are  assumed  to  lie.  Two  vectors  that  define  the  plane,  PI  P2  and  P2  P3, 
are  then: 

PfP2  = <X,a  - XM)i  + <Y„  - YP1)i  + (Z,a  - Z„)k  (113) 

P2P3  = <X,S  - X,a)i  + (Y„  - Y„)j  + (Z„  - Z,a)k  ( * ' 


The  normal  to  the  plane  defined  by  these  vectors  is  given  by  the  cross 
product: 

N = PI  P2  x P2  P3  (115) 

which  in  matrix  notation  is: 

k 
I 

N = 


i J 

(X„  " XM)  (\^a  - Y^)  (Zpa  - ZPj) 


l?a  * n' 


(X„  - X,a)  (Y„  - Y,a)  (Z„  - Z,a) 


M 1 P3' 


'»s 


(116) 


or: 


N = 


i 

Xfl 


k 

ai  ^ai 


(117) 


1 sa  J 38 


where  X,,,  Y^  . , . are  the  corresponding  quantities  in  Equation  116.  Ex- 
panding the  determinant  in  Equation  117  gives: 

N - (Ygj Zj,  - Y^ZgjJi 

' tt^Z*,-  X^Z,,))  (118) 

+ «8iYm-  X^Y^k 

The  unit  vector  normal  to  the  plane  will  be: 

NU  = NUX  1 + NU  Y j + NUZ  k (1 19) 

where: 

NUX  * <YaiZM-  Y^Z^J/NU 
NUY  * (XjjZ*  - Xj.Z^J/NU 
NUZ  = (X*  - X^  Yal)/NU 

NUa  = (YaiZ„  - YMZa)a  + (Xa,Z„  - XMZaiP  + (X3l  Y„  - X,aYai)1 * 3 
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NUMERICAL  INTEGRATION  OF  FLUX  DENSITY 


Th*1  total  magnetic  flux  normal  to  the  loop  is  determined  by  a numerical 
integration  process.  The  process  is  shown  in  Figure  42.  The  plane  is  divided 
vertically  and  horizontally  into  12  equally  spaced  strips.  For  this  discussion, 
vertical  will  mean  the  direction  of  point  1 to  point  2 or  point  4 to  point  3,  and 
horizontal  will  mean  the  direction  of  point  PI  to  point  P4  or  point  P2  to  point  P3. 


Figure  42.  Integration  Techniques  for  Flux  in  Plane 


The  12  strips  define  13  lines  vertically,  the  intersections  of  which  define 
169  points  (13  y 13),  of  which  point  PL  (9,  5)  is  shown.  The  flux  density  at 
each  point  is  evaluated  and  then  integrated  numerically  along  each  of  the  ver- 
tical strips.  The  resultant  13  values  are  integrated  horizontally  to  obtain  the 
total  magnetic  flux  linking  the  plane.  The  integration  process  used  is  called 
Weddle's  rule  (Ref.  11)  and  is  based  on  fitting  a sixth  order  polynomial  to  the 
array  of  points  and  then  integrating  the  resultant  polynomial.  The  result  is: 

F H • AX  = Jq  AX  (11,  +•  Slip  + 11,  + 6II4  + 11,  + 511.  ♦ IIJ  (120) 

In  program  APERTURE,  two  such  polynomials  are  fitted  to  the  13  points, 
giving  the  following  coefficients: 

1,  5,  1,  0.  1.  9,  2.  5,  1.  6.  1,  5.  1 
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DEVELOPMENT  Of  COMPUTER  PROGRAM 


PROGRAM  DESCRIPTION 

A program  listing  for  APERTURE  is  given  in  Figure  43.  The  listing  shown 
consists  of  the  MAIN  portion  of  the  program  and  two  subroutines:  SHAPEFAC, 
which  is  used  once  during  the  running  of  MAIN;  and  MAGFLD,  which  is  used 
repetitively  during  the  running  of  MAIN.  Figure  44  is  an  elementary  flow 
chart  for  MAIN. 

Hie  program  first  reads  the  input  data  from  a file,  the  name  of  which 
will  be  requested  during  the  program.  It  next  determines  the  effective  dipole 
moments  presented  by  the  aperture,  in  both  the  X and  Y directions.  X and 
Y are  taken  to  be  oriented  along  the  major  and  minor  axes  respectively  of 
the  aperture. 

Hie  program  next  tabulates  the  magnetic  field  intensity  at  the  desired 
points  of  the  region  beyond  the  aperture.  This  tabulation  may  or  may  not 
include  the  effects  of  a reflecting  surface  behind  the  aperture.  This  portion 
of  MAIN  uses  the  subroutine  MAGFLD  to  calculate  the  field  intensities  at  the 
point  under  consideration.  Should  the  tabulation  of  field  intensitiec  not  be 
desired,  this  portion  of  the  program  is  bypassed. 

The  program  then  goes  on  toward  the  calculation  of  the  flux  that  passes 
through  a four-sided  loop.  Hie  loop  is  defined  by  theX,  Y,  Z coordinates 
of  the  four  points  making  its  comers.  The  first  three  points  are  used  to 
define  the  plane  of  the  loop;  the  fourth  point  is  assumed  to  lie  in  this  same 
plane.  After  reading  the  coordinates  of  the  defining  points,  the  program  cal- 
culates the  field  intensity  at  169  points  over  the  surface  of  that  loop.  A nu- 
merical integration  of  the  field  intensity  at  these  169  points  is  then  performed 
to  obtain  the  total  field  intensity  and  total  flux  passing  through  the  loop.  After 
completing  the  calculation  of  total  flux  through  the  first  loop,  the  program 
reads  the  coordinates  of  additional  points,  and  calculates  the  flux  through 
such  planes  as  may  be  defined.  Hie  program  continues  to  run  in  this  manner 
until  no  further  loops  are  encountered.  If  desirable,  this  portion  of  the  pro- 
gram also  may  be  bypassed. 

Input  data  --  long  form  — for  the  APERTURE  program  are  shown  in 
Figure 45;  Figure  46  is  the  short  form. 

MAIN  Program 

Before  starting  the  detailed  description  of  the  MAIN  program,  the  user 
should  refer  to  Figure  47,  which  gives  the  terminology  by  which  the  aperture, 
the  external  magnetic  field,  and  the  point  under  consideration  are  described. 

Hie  X,  Y,  Z coordinates  of  both  the  aperture  and  the  point  under  consideration 
are  given  with  respect  to  a reference  set  of  axes.  The  plane  containing  the 

*This  listing  is  for  a program  that  will  be  run  on  the  General  Electric  Time 
Sharing  computer.  A program  listing  for  the  CDC6600  machine  is  included 
in  Appendix  III,  "Program  Listings  for  CDC6600  Computer.  " 
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APERTURE  I A*  37E5T  00/05/75 

IOOOC  APERTURE A PRC  SRAM  THAT  CALCULATES  THE  MAGNETIC  FIELD  THAT 

tOIOC  PASSES  THROUGH  AN  APERTURE*  FA  FISHER  BLDG  9*209 

IOOOC  GENERAL  ELECTRIC  COMPANY  tOO  WOOOLAWN  AVE  PITTSFIELD* MASS  01201 

I030C  PHONE  (41 3) -494-4300 

1040C  DEVELOPED  UNDER  CONTRACT  F3361 1 -74-C-3068  USAF  FLIGHT  DYNAMICS  LB 
IOSOC  THE  PROGRAM  READS  DATA  FROM  AN  EXTERNAL  FILE*  THE  NAME  OF  MUCH 
1060C  WILL  BE  REOUESTED  DURING  EXECUTION.  THE  INPUT  OATA  FILE  SHOULD 
1070C  BE  CONSTRUCTED  AS  FOLLOWS! 

IOOOC 
I090C 
ItOOC 
II  IOC 

neoc 

II30C 
II40C 
II50C 
!1«0C 
II  TOC 

tteoc 

II90C 
IOOOC 
IP  toe 
ioo^c 
1C30C 
I040C 
1050C 
IOSOC 

1070C  LI  AND  LO  ARE  THE  LENGTHS  IN  METERS  OF  THE  AXES  OF  THE  ELLIPTICAL 
IOOOC  APERTURE*  L I "MAJOR  AXIS  AND  » LO-MINOR  AXIS* 

ie90C  AN  AH  IS  THE  ANGLE  THAT  THE  MAJOR  AXIS  OF  THE  APERTURE  MAKES  WITH 
1 30 OC  THE  X AXIS*  0 DEGREES  IS  PARALLEL  THE  THE  POSITIVE  X AXIS* 

13I0C 

iseoc  NEXT  IS  THE  strength  in  amperes  per  meter  of  the  external  field 

I330C 

I340C  ANGH  WITH  RESPECT  TO  THE  X AXIS*  0 DEGREES" PARALLEL  TO  X-AXIS* 
I3S0C  Dt"l *YE5-THERE  IS  A REFLECTING  SURFACE  PARALLEL  TO  THE  APERTURE* 
I3S0C  D1«0"N0  REFLECTING  SURFACE* 

I370C 

I300C  D2"Z  COORDINATE  OF  THE  REFLECTING  SURFACE*  ENTER  DUMMY  VALUE  IF 
I390C  Dt"0* 

I400C 

1 41 OC  D3» I "Y ES-CALCULATE  THE  FIELDS  OVER  A PRESCRIBED  VOLUME  INSIDE* 
I400C  D3"0»N0-SKIP  THIS  CALCULATION. 

I430C 

I440C  ZPA«Z  COORDINATE  AT  WHICH  CALCULATION  SHOULD  START 

Figure  43.  APERTURE  Program  Fisting  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  1 of  9) 


LINE  NUMBER  10  XA*YA*ZA 

00  LI*LO* ANAH 
30  NEXT*  ANGH 
40  Dl  * DO 
SO  D3 

60  ZPA*ZPB*ZPC 
70  YPA*YPB*YPC 
80  XPA*XPB* XPC 
90  D4 
100  DS 

no  Pxi*PYt«pzi*pxe*PYO*pze 
too  PX3*PY 3*PZ3*PX 4* PY4»PZ4 

(LINE  NUMBERS  NEED  NOT  BE  IDENTICAL  TO  THOSE  ABOVE) 


XA*YA*ZA  ARE  THE  COORDINATES  IN  METERS  OF  THE  CENTER  OF  THE 
APERTURE*  IT  IS  LOCATED  IN  A PLAN*  PARALLAL  TO  THE  XY  PLANE 
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I4S0C  ZPB*Z  COORDINATE  AT  WHICH  CALCULATION  SHOULD  END 
I460C  ZPC*Z  INCREMENT  SIZE 

I470C  YPA# YPB# YPC#XPA#XPB#XPC  ARE  SIMILAR  F#R  X AND  Y COORDINATES 

! 48 OC  ENTER  DUMMY  VALUES  IF  D3-0 

I490C 

1SOOC  D4* I "TABULATE  FIELD  IN  SPHERICAL  COORDINATES* 

ISIOC  D4-0-TABULATE  IN  RECTANGULAR  COORDINATES* 

ispoc 

1530C  D5* I "YES-CALCULATE  THE  FLUX  LINKING  A LOOP 
IS40C  D5-0-N0-SKIP  THIS  CALCULATION. 

IS50C 

1S60C  PX1.PY1# PY4.PZ4  ARE  THE  COORDINATES  OF  FOUR  POINTS  THAT 

1S70C  DEFINE  THE  LOOP*  THEY  MUST  00  AROUND  THE  LOOP  IN  CONSUCUTIVE 
I WOC  ORDER*  ADDITIONAL  LOOPS  MAY  BE  DEFINED  BY  ADDITIONAL  DATA  IN 
1590C  THE  SAME  FORMAT.  DUMMY  VALUES  ARE  NOT  REOUIRED  IF  D5*0 

1600C  

1610  FILENAME  INFILE 

1620  REAL  LI#L2#NU1#NU2#NU3*NU#NUX#NUY#NUZ 
1630  DIMENSION  HN! 1 3# 1 3) 

1640  DIMENSION  T86A(I3) 

I6S0  DIMENSION  PATHA<I3> 

1660  10  PRINT  15 

(670  15  FORMAT!”  INPUT  FILE  NAME**) 

1680  20  INPUT# INFILT 
1690  30  PRINT  I IS 

I700C  CARRIAGE  CONTROL  FORMAT  STATEMENTS 


FORMAT! 1H-) 
FORMAT! I HO) 
FORMAT! IH  ) 
FORMAT! 1H 6) 
FORMAT! IH*) 


1710  110 
1720  115 
1730  120 
1740  122 
1750  123 
I760C  OUTPUT  DATA  FORMATS 
1770  130  F0RMAT!6EI2*3) 
1780C  DATA  HEADING  FORMATS 
1790  140  FORMAT!”  APERTURE 
FORMAT!” 

FORMAT!” 

FORMAT!”  APERTURE 
FORMAT!” 

FORMAT!”  APERTURE 
FORMAT!”  EXTERNAL 


1800 
1810 
1820 
1830 
18  40 
1850 
I8  60 
1870 
I860 
1890*  ” 
1900  188 
1910  190 
1920  192 


145 

ISO 

155 

160 

165 

170 

175 

180 

185 


COORDINATES— X«”#IEie. 3#”  METERS”) 
Y*”#IEIB*3#”  METERS”) 
Z*”#  IEt£*3»”  METERS”) 
DIMENSIONS- -MAJOR  AXIS"”#  IEIt*3#” 
MINOR  AXIS"”#  1E12.3#” 
INCLINED”# I EtO*3#”  DEGREES  FROM  X 
MAGNETIC  FIELD*”# I El «• 3#”  AMPERES 


METERL”) 

METERS”) 

AXIS”) 

PER  METERP 
AXIS”) 


FORMAT!”  AND  INCLINED”# IEIB*3#”  DEGREES  FROM  THE  X 
FORMAT!”  THERE  IS  NO  REFLECTING  SURFACE**) 

FORMAT!”  THERE  IS  A REFLECTING  SURFACE  LOCATED  AT  Z*”#IEIB*3# 
METERS”) 

FORMAT!”  LOOP  NUMBER  ”. I 5) 

FORMAT!”  LOOP  AREA*”# I El 0*3#”  SGUARE  METERS”) 

FORMAT!”  TOTAL  FLUX*”#  I El f. 3.  " WEBERS”) 


Figure  43.  APERTURE  Program  T ?sting  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  2 of  9) 
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If  5 
200 
210 
220 


PRINT 

PRINT 

PRINT 

PRINT 

PRINT 


1930 
1940 
19  50 
I960 
1970 
1980 
1990 
2000 
2010 
2020 
2030 
2040 
2050 
2060 
2070 
2080 
2090 
2100 
2110 
2120 
2130 
2140 
2150 
2160 
2170 
2180 
2190 
2200 
2210 
2220 
2230 
2240 
2250 
2260 
2270 
2280 
2290* 

2300  68T8 
2310  2206 
2320  2207 
2330* 

2340  2208 


FORMATC"  8UT  8P  DATA") 
FORMATC V) 

F8RMAT<"  POINT  X 
FORMATC 1 5#  3EI 2. 3) 


READC  INF  1LE#  200#  END* I96Q)LINE»XA#  Y A#  2 A 
READ  C 1NF1LE# 200* END*1960)L1NE#L l#L2# ANGA 
PRINT  I 40# X A 
1 45*  YA 
I 50»  ZA 
I 55, LI 
160.L2 
1 65*  AN6A 

READC  1NF1LE#  200*  END-  I960>L1NE»H£XT»  ANON 

PRINT  115 

PRINT  1 70# HEX T 

PRINT  I 75,  ANGH 

PRINT  115 

READC  1NF1LE# 200#  END"  1 9 60) LINE,  Dl  # D2 

275  1FC 01)280# 280# 290 

280  PRINT  180 

PRINT  115 

28  5 68T8295 

290  PRINT  185# 02 

PRINT  IIS 

295  CONTINUE 

READC INFILE# 200# END- 19 60)LINE# 03 
READC 1 NFILE, 200, END* I 9 60)L I NE# ZPA# ZPB* ZPC 
READC  INFILE,  200#  END*  I960)L1NE#  YP4#  YPB#  YPC 
READC 1NF1LE# 200# END-I960)LINE»XPA#XPB#XPC 
READC  1NF1LE* 200,  END- 19  60)L  1 NE,  D4 
PI-3.141  59  2 6 5 

CALL  SHAPEFACCL1#L2#A1 l#A22) 

1FCD3)  19  SO#  1950.2200 

2200  1FC 04) 220 1 #220 1 #2206 

2201  PRINT  2202 

2202 


2350 

2360 

2370 

2380 

2390 

2400 


FORMATC"  X 

H-Y  H-Z") 

2204 

PRINT  2207 
F8RHATC " X 

LAT  LONG") 

PRINT  120 
GOTO  2450 
2204  PRINT  120 
2450  CONTINUE 
J I -I F1XC CZPB'ZPA) /ZPC)  ♦ I 
J2-I FIXC CYPB"YPA)/YPC>* 1 
J 3- 1 FIXC  CXPB-XPA) /XPC) ♦ I 


Z"> 


H-X 


HTOT 


Figure  4 3. 


AI’KKTUKPJ  Program  Listing  for  the  General  Eleetric 
Time  Sharing  Computer  (Sheet  3 of  0) 
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2410  Dt  If  SO  1 3a1*  J3«  1 
2420  DOlf SO  I2a1* J2* 1 
2410  DOlf SO  lla1*JI*| 

2440  XPt«XPA4(13-1>*XPC 
2450  fPI»YPA*<!2-1>*YPC 
2460  ZPt>ZPA4(11-1>*ZPC 
2470  1002  CONTINUE 
2480  tiff  CONTINUE 

24f0  1 750  CALL  NAGFLD<ANGA*ANGH*XPl*YPl*ZPl*XA«YA*ZA#HEXT#A1  1#A22# 
2500*  HPX 1 # HPY I *HPZ1 *L1 *L2* D1*  D2> 

2510  3370  I F< D4) 1210# 1210#  4000 

2520  1210  PRINT  1220*XP1#VPI#ZP1*HPXI*HPYt*HPZt 

2530  GOTO  If  SO 

2540  1220  F0RMAT(6E12*3) 

2550  4000  D*SORT < HPX I *HPX 1 *HPZ I *HPZ 1 > 

2500  4002  I F(  ABS<  HPY I > "ABS<  D)  ) 400  4#  400  4*  4010 
2570  4004  AN61«f0-S7.2f577fS*<ATAN<HPYI/D>> 

2560  4006  «T0  4012 

2Sf0  4010  AN6 1 * 57 • 2f  577f  5f<  ATAN<  D/HPY I ) > 

2600  4012  1 F< ABSCHPX 1 ) -ABS(HPZ I ) ) 40 1 4# 40 1 4# 4020 
2610  4014  ANG2afO~S7*2f S77f 5*(ATAN<HPX1/HPZ 1 >> 

2620  4016  GOTO  4030 

2630  4020  ANG2* 5T*2f 577f 5*< ATANCHPZ 1/HPX1 > > 

2640  4030  IP(HPYI)  4050#  4050*  4040 

2650  4040  GOTO  4110 

2660  4050  ANG1 a180~ANG1 

2670  4110  CONTINUE 

2680  4120  I F<HPX  1)4180#  41 30#  41 30 

26f0  4130  I F<HPZ  1)41 60# 41  40#  41  40 

2700  4140  ANG2«ANG2 

2710  4150  GOTO  4215 

2720  4160  ANG2«-ANG2 

2730  4170  GOTO  4215 

2740  4180  I PCHPZ 1 ) 421 0#  41f 0#  41f 0 

2750  41f 0 ANG2a180*ANG2 

2760  4200  GOTO  4215 

2770  4210  ANG2a-(l80>ANG2) 

2780  4215  CONTINUE 

27f 0 HPTaSORT<HPX I 4HPX I ♦HPY 1 *HPY I 4HPZ 16HPZ I > 

2800  4230  PRINT  1220*XP1 #YP1# ZP1#HPT# ANG1 # ANG2 
2810  If SO  CONTINUE 
0820  PRINT  110 

2830  REAO(  INFILE#  200#  END* If  60) LINE*  D5 
2840  2100  I F< 05)1400*1400*  If  55 
8850  IfSS  CONTINUE 
2860  JX*0 

2870  If 57  CONTINUE 

2880  READt INFILE* 200*  ENDalf 60)LINE#PX!#PY  t#PZt#PX2«  PY2#PZ2 

Figure  43.  APERTURE  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  4 of  9) 


88*0  READ<INP1LE#800#END«1*60>LINE«PX3#PY3#PZ3#PX4#PY4*PZ4 

8*00  C THESE  ARE  THE  SIDES  «P  THE  QUADRILATERAL 

8*10  8110  CONTINUE 

8*80  X81*PX8-PX1 

8*30  X38«PX3-PX8 

8*40  X43-PX4-PX3 

8*50  X14-PX1-PX4 

8*60  Y81-PY8-PYI 

8*70  Y38*PY3-PY8 

8*80  Y 43*PY  4*PY  3 

8**0  Y14-PY1-PY4 

3000  Z81«PZ8*PZ1 

3010  Z38-PZ3-PZ8 

3080  Z43«PZ4-PZ3 

3030  Z14»PZ1-PZ4 

3040C  THIS  IS  A DIAGONAL  OP  THE  QUADRILATERAL 
3050  X31«PX3*PX 1 
3060  Y31*PY3“PY1 
3070  Z31*PZ3-PZ1 

3080  T81«SQRT(Xei*xei«Y81*Yei4Z8t*Z8l> 

3090  T38*SQRT<  X 38*X 38* Y 3fi*Y 38+Z  38*Z  38) 

3100  T43*SQRT(X  43*X 43*Y 43*Y  43+Z  434Z43) 

3110  T1 4*S0RT<X 1 4*X 1 4+Y 1 44Y 1 4+Zl  4*Z1  4) 

3180  T31«S0RT(X31«X31*Y31*Y31*Z31*Z3t> 

3130  SlXT81«T38*T31>/8 

3140  A1*SQRT< Sl*< Sl-T81>*<  SI*T38>*< S1"T31 ) ) 

3150  S8«(T43*T14»T31)/8 

31 60  A8>S0RT( S8*C S8-T43>*< S8-T1 4>*< S8-T3I ) > 

3170  AREA-A1+A8 

3180C  THESE  ARE  THE  MIDPOINTS  OP  THE  ENDS  OP  THE  QUADRILATERAL 

31*0  XPMl«PXl*X8t/8 

3800  YPM1-PY1 >Y81/8 

3810  ZPM1«PZ 1*18 1/8 

3880  XPM8>PX4*X43/8 

3830  YPM8>PY4«Y43/8 

3840  ZPM8*PZ  4-Z  43/8 

3850  XPM81>XPM8-XPM1 

3860  YPM6 1 *YPM8* YPM 1 

3870  ZPM8 1 *ZPM8" ZPM 1 

3880  TPM« SO  RT( X PM8 1 *X  PM8 1 ♦ YPM8 1*  YPM8 1 *ZPM8 1 *ZPM8 1 > 

38* OC  THESE  ARE  THE  COMPONENTS  OP  THE  NORMAL  VECTOR 

3300  NUl«Y8l*Z38-Y38*Z81 

3310  NU8«*X81*Z38*X38*Z81 

3380  NU3«X81*Y38-X38*Y81 

3330  NU*3QRT(NUl*NUl+NUt6NU8+NU3*NU3) 

3340C  THESE  ARE  THE  COMP ONON ENTS  OP  THE  UNIT  NORMAL  VECTOR 
3350  NUX*NU1/NU 
3360  NUV*NU8/NU 

Figure  43.  APERTURE  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  5 of  9) 
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3370  NUZ=NU3/NU 
3390  3550  CONTINUE 
3390  3560  DO  3880  N2a|»13»t 
3400  3570  DO 3880  Nl«l«13»l 
3410  XP5aPX!*Xei*<Nl-l>/ie 
3480  YP5aPYl*Y21*<Nl-l>/l2 
3430  ZP5»PZl*Z8l*<Nl-l)/ie 
3440  xpr»apxe+x32*<Ne-i>/ie 
3450  YP6-PY2*Y32*<N2-1>/12 
3460  ZP6»PZ2«Z32*<N2-1>/12 
3470  XP7-PX4-X43*<N1-1>/12 
3480  YP7«PY4-Y43*<N1-1>/12 
3490  ZP7-PZ4-Z434CN1-D/12 
3500  XP8»PXI-X14*<N2-1>/12 
3510  YP8-PY1-YI 44<N£-1>/12 
3520  ZP8»PZ1-Z14*<N2-1>/12 
3530  X75-XP7-XP5 
3540  Y75-YP7-YP5 
3550  Z75«ZP7-ZP5 
3560  X86"XP8-XP6 
3570  Y86-YP8-YP6 
3580  Z86-ZP8-ZP6 

3590  T8 6* SORTC X8  64X8  6« Y8  6*Y8  6+Z8  6*ZB 6) 

3600  XPL«XP8-X864<NI-1>/12 
3610  YPL-YP8-Y864CN1-D/12 
3620  ZPL»ZP8-Z864(N|-1>/12 

3630  CALL  MAGFL D( AN  GA» ANGH#  XPL# YPL#  ZPL*  X A* YA# ZA»  HEX  T»  A1 1 t A£2# HPX I 
36406  * HPY I *HPZ 1 »L  I #L2#  Dl»  D2> 

3650  HNP-HPX  1 4NUX*HPY 1 4NUY4HPZ 1 4NUZ 

3660 C THESE  ARE  THE  HN'S  AT  THE  VARIOUS  PRINTS  RP  THE  QUADRILATERAL 
3670  HNCN1#N2)»HNP 

3680C  THESE  ARE  THE  DISTANCES  TRP  TR  DRTTRN  ALONG  THE  QUADRILATERAL 

3690  T86A(N2)«T86 

3700  3880  CONTINUE 

3710  3890  0039 50  N2-l»13»l 

3780  DELTA 1*T86ACN£>/ 12 

3730C  THIS  EVALUATES  FLUX  ALONG  THE  LINES  IN  THE  DIRECTION  P1**»P£ 
3740C  AND  P4--*P3 

3750  PATH»HN<  1 » N2)  ♦HNC  8#  N2)  6 54HNC  3*  N8 ) «HN<  At  N8)  6 6>HN(  5#  N2)  4HNC  6*  N8» * 5* 
37606  HN(7#N2)62»HNC8aN2)*5+HN<9#N2)4KNC 10*N2>464 
37706  HNC 1 1»N£)4HN( 12»N£>*54HN( I3*N2) 

3700  PATHA(N£)«0«3*DELTA1*PATH 
379:  39  50  CONTINUE 
3800  DBLTA2«TPM/12 

38 IOC  THIS  EVALUATES  THE  RESULTANT  FLUK  IN  THE  DIRECTION  P1-->P4 
3880C  AND  P8-->P3 

38  30  HTOT-PATHAC 1 > 4PATHAC £>  4 5+PATHAt 3) ♦PATHAC 4> * 6*PATHA( 5) ♦ 

38406  PATHAC 6>454PATHA(7)4£«PATHA<8)*54PATHA(9>4PATHA( 10)464 

Figure  43.  APERTURE  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  f of  9) 


38  50*  PATHA!  1 1 ) ♦PATHA!  18)*5«PATHA!  13) 

3860  HT0T«0  « 3* DEL TA8*H TiT 

3870  BT0T*HT0T*4*PI*lE-7 

3880  JX«JX«1 

3890  4315  PRINT  188* JX 

3900  PRINT  180 

3904  PRINT  810 

3906  PRINT  180 

3910  J»1 

3980  PRINT  820* J*PX1*PY1«PZI 
3930  J*2 

3940  PRINT  880* J* PX8* PY8* PZ8 
3950  J-3 

3960  PRINT  880* J*PX3*PY3*PZ3 

3970  J>4  i 

3980  PRINT  820* J*PX4*PY4*PZ4 

3990  PRINT  180 

4000  4380  PRINT  190*AREA 

4010  4330  PRINT  198*BT8T 

4080  4340  PRINT  110 

4030  68T8  1957 

4040  4360  PRINT  110 

4050  1960  PRINT  195 

4060  1400  STOP 

4070  END 

4080  SUBROUTINE  SHAPEFAC(L1*L2*A1 1.A28) 

4090  REAL  L1*L8 
4100  PI *3. 141  59865 
4110  EI*1-CL8/L1)**8 
4180  E8-SQRTCE1) 

4130  lF!L2/Li«0)G0T0  3130 

4140  IP(L8/L1»1)  GOT*  3160 

4150C  00000991 09V0000000000000C 000000 

4160C  CELI(  1*  £2)  AND  CELI(8*E8>  ARE  MATH  LIBRARY  ROUTINES  THAT 
41  TOC  EVALUATE  THE  ELLIPTIC  INTEGRALS  OP  THE  FIRST  AND  SECOND  KINDS. 

4180  Y1»CEL1!1*E2) 

4190  Y2*CELI(8*E8) 

480 OC  0000000000000090900909990000000 

4810  C0Nl>8*PI*(Ll/2)**3/3 
4880  A11«C0N1*E1/(Y1*Y8) 

4830  A88«C0N1*E1*( I”E1 )/!Y2”! 1”EI >*Y 1 > 

4840  A33«C0N1*(1-EI)/Y8 
4850  RETURN 

1 1 1ft  PRINT  11  iA  X 

4870  3140  FORMAT! " L8/LI  IS  NEGATIVE*  THIS  IS  AN  ERROR") 

4880  RETURN  I 

4890  3160  PRINT  3170 

4300  3170  FORMAT!"  L8  IS  LARGER  THAN  LI*  THIS  IS  AN  ERROR") 

Figure  43.  APERTURE  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  7 of  9) 
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4310  STOP) END 

4320  SUBROUTINE  HAGFLDC AN6A*  AN6H/XP1/ YP1*ZPI*XA*YA*ZA*HEXT*  A1 1* 

43304  A22/HPX1/HPYI/HPZ1/L1/L2/D1/D2) 

4340  REAL  KI*K2*KS*K6*L1*L2 
4350  N»0 

4360  RAD»»57. 2957795 
4370  Pla3>  1 41 S926S 
4380  K I "COSC ANQA/RAO) 

4390  K2-SINCAN6A/RAD) 

4400 C CALCULATIBN  OF  SHIFTED  COORDINATES  OF  POINT  UNDER  INVESTIGATION 

4410  XP2-XP1*K1*YPI*K2 

4420  VPB*-XPI*K2*YP1*K1 

4430  ZP2-ZPI 

4440  HPX2*0 

4450  HPY2*0 

4460  HPZ2*0 

4470  ZAA«ZA 

448 OC  CALCULATION  OF  DISTANCES  FROM  APERTURE  TO  POINT  UNDER  STUDY 
4490  9130  XC*XP2-XA 
4500  9140  YC«YP2-YA 
4510  9150  ZC*ZP2-ZAA 

4520  9160  Cl>XC*XC*YC*YC*ZC*ZC*LI«Ll/4 
4530  C2*XC«XC*YC«YC*ZC*ZC*L2*L£/4 

4S40C  CALCULATION  OF  H FIELD  PARALLEL  TO  AXES  OF  APERTURE 
4550  CONl*l • 43833 
4560  C0N2» 1*804688 
4570  C0N3*2.094727 

458 OC  CALCULATION  OF  FIELD  PARALLEL  TO  AXES  OF  APERTURE 
4590  ANAH>(AN0H-AN6A)/RAD 
4600  HHAJ*H  EX  T*COS( ANAH > 

4610  HMIN«HEXT*S1N<ANAH> 

4620 C CALCINATION  OF  ROTATED  COMPONENTS  OF  MAGNETIC  FIELD 
4630  KS«-A11*HNAJ/(4*PI*L1> 

4640  K6»A22«HMIN/<4»PI*L2> 

4650  C3«*«*LI 

4660  91  65  IF(ABS( C3>*1  E*S>  9100*9170*9170 

4670  9170  F3*C3*C3 

4680  F7"1»£*F3*59F3*F3*7*F3*F3*F3 

4690  GOTO  9190 

4700  9180  F7«l 

4710  9190  C5«C3/C1 

4720  9195  k F<  ABS< CS>- 1 E” 5)  9210*9200/9200 
4730  9200  C7«C5*C5 

4740  Fl«C5*t  1*C0N1*C7«-C0N2*C7*C7*C6N3*C7*C7*C7) 

4750  GOTO  9220 

4760  9210  F1»C5 

4770  9220  C4»YC*L2 

4780  9 230  IF'ABSCC4>-IE»5>  9 2 50*9  2 40/92  40 

Figure  43.  APERTURE  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  CSheet  8 of  9) 
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4790  9240  rS»C4*C4 
4boo  f8*i»2*f5*5*f5*f50*F5*fs*F3 
4810  6# TV  9fi<0 
4820  92S0  FS*1 

JK1A  ffliA  riaCi/CP 

4840  9 270  IF<ABSCC6>-IE-S>  9300*9280*9280 
4850  9280  C8«C6*C6 

4860  F2«C6*<1*CVN1*C8*C8N2*C8*C8*C8N3*C8*C8*C8> 

4870  08 TV  9320 
4880  9300  F2-C6 
4900  9320  CONTINUE 
4910  8l«3*K5*XC/CI**l • 5 
4920  82*39K6*XC/C2**I • 5 
4930  63»KS*L|/C!**1«S 
4940  G4»3*KS*YC/Cl*4 1 • 5 
4950  Q5*3*K6*YC/C29*1 • 5 
4960  G6>K6»L2/C2**1«5 
4970  G7»3*KS*ZC/C1**I>5 
4980  68»3*K6*ZC/C2**1«S 
4990  HPX»6l*Fl*G2*F2,,G3*F7 
5000  HPY«64*F1485*F2-86*F8 
5010  HPZ«67*Fl»68*F2 
5020  9680  I F( 01  > 9760*9760*9682 
5030  9682  HPX2«HPX2*HPX 
5040  HPY2»HPY2*HPY 
*050  HPZ2«HPZ2*HPZ 

5060  9 690  IFCABSCHPX)a«00S6ABS(HPX2>)9700»9  700*9720 
5070  9700  I FCA8S(HPY)~ *00 St'ABSCHPYB)  >9710*9710*9720 
5080  9710  1F<ABSCHPZ) ~*00$4ABS(HPZ2> >9770*9770*9720 
5070  9720  I FCN- 10)9730*9730*9770 
5100  9730  N»N*l 

5110  9740  ZAA«ZAA~< C •! >**N>*2*N*D2 

5120  6VTV9I30 

3130  9760  HPX2*HPX2*HPX 

5140  HPY2-HPY2«HPY 

3150  HPZ2»HPZ2*HPZ 

5155  9 770  CONTINUE 

5160C  CALCULATION  OF  CVHPONENTS  O'  NAVNETIC  FIELD  ROTATED  BAC* 

51  TOC  TV  THE  REFERENCE  AXE5 
5180  HPX 1 »MPX2*K 1 "HPY24K2 
5190  HPY 1 •HPX2*K2^HPY24K I 
5200  HPZ1*HPZ2 
5210  RETURN 
5220  5230  STOP) END 

Figure  43.  APERTURE  Program  Listing  for  the  General  Electric 
Time  Sharing  Computer  (Sheet  9 of  9) 
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ft 


No 

(D3  = 0) 


Lilies  1970 


Figure  44.  Main  Program  --  Elementary  Flowchart 


Coordtoitu  of  the  Aperture 


tt  "jr-  ryratjrz-  -\~FW  ~ . <*yv  «»m^  '»r  k "7*^  ’» 


- r^T'fTr^  ,-*?T 


in-  X Coordinate.  Y Coordinate  Z Coordinate 

(XA)  (YA)  (ZA) 

PiwwrtOM  tod  Orientation  of  the  Agtrtw 

Orientation  of  Major  Axis 

Major  Axis  Length  Minor  Axis  Ufi|th  with  Respect  to  the  X-Aais 

20b , , (ANAH) 

(Li)  (Lit  (ANAH) 

Magnitude  end  Direct  too  o:  the  Externxl  Magnetic  Field 

Field  Strength , Orientation  with  Reaped  to  X-Axis 
JOb  ~ (HE*TT  fANGK)  


40b 


Reflecting  Surface 


le  there  a reflecting  surface?  For 
Yaa  antar  1 (one)  and  for  No  enter 
0 (zero) 


(DO 


How  far  from  the  aperture  la  the 
reflecting  surface?  if  there  ia  no 
reflecting  surface  enter  the  dummy 
number  1000 
(D2) 


SOb 


50b 


70b 


eob 


»0i 


Tabulat  ion  of  Field  Strength 

Do  you  want  a table  of  the  ftetd  atrengtha 
over  the  interior  regioof  For  Yes  enter 

1 (ooe)  and  for  Wo  enter  0 (zero) 

(D3) 

Volume  Orer  Which  Fielda  -xe  to  be  Calculated 


Z Dimension 

s)  Start  st,  b)  End  at . clAnd  increment  in  Steps  of  (ZPC) 

(ZPA)  (ZPB)  (ZPC) 

Y Dimension 

a)  Start  st  . b)  End  st  . c)And  Increment  in  Steps  of , 

(YPA)  fYPB)  fYPC) 

X Dimension 


a)  Start  st . 
(XPA) 


b)  End  st  , c)And  Increment  in  Steps  of . 

(xpb)  ixpcS 


[On  lines  60,  70,  and  00  enter  the  dummy] 
numbers  0 (xero)  in  each  of  the  locations  I 
if  you  don’t  want  the  fields  tabulated  J 


Oitput  Format 

Do  you  want  the  tabulation  in  rectangular 
or  spherical  coordinates?  Enter  1 (one) 
for  spherical  or  0 (zero)  fhr  rectangular. 
Enter  a dummy  rtloe  if  no  tabulation  la 
desired. 

fin) 


Total  Flux  In  a Loop 

Do  you  want  to  determine  the  total  flux 
linking  a loop?  Enter  l (one)  for  Yea  and 
0 (aero)  for  No. 

(D5j 


1 i 0b 


120b 


Enter  the  X,  Y , Z coordinates  of  tbs 
four  points  defining  the  loop.  They 
should  go  in  sequence  around  the  loop. 
All  tour  points  must  lie  In  the  same 
plane.  If  a loop  calculation  la  not  de- 
sired the  following  linea  may  be  left 
blank,  in  which  case  the  program  makes 
an  automatic  atop. 


Point  1 

X , Y , Z 
(PX1)  (PYi)  (PZl) 

Point  3 

C'XJ)  TpvT)  TpzTi' 


Point  2 

X . Y , Z . 
(PX2)  TtH)  (PZ2> 

Point  4 

(PX4t*  W\d'  iPU) ' 


Successive  lines  following  the  format 
of  110  and  120  may  be  used  to  enter  the 
defining  points  for  other  loops  The  pro- 
gram will  automatically  stop  when  It  runs 
out  of  data. 

Figure  45.  Input  Data  for  Program  APERATURE 


Long  Form 


10b 


20b 

(XA) 

» 

(YA) 

fLl) 

(L2) 

30b 

» 

(HEXT) 

(ANGH) 

40b 

* 

(Dl) 

(D2) 

50b 

(D3) 

60b 

» 

i 

(ZPA) 

(ZPB) 

70b 

» 

» 

(YPA) 

(YPB) 

80b 

(XPA) 

(XPB) 

90b 

(D4) 

100b 

(D5) 

110b 

* 

(PX1) 

(PY1) 

120b 

1 

J 

(PX3) 

(PY3) 

(ZA) 


(ANAH) 


(ZPC) 


(YPC) 


(XPC) 


(PZ1)  (PX2)  (PY2)  (PZ2) 

(PZ3)  (PX4)  7PY4)  (PZ4) 


Figure  46.  Input  Data  for  Program  APERATURE  --  Short  Form 
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Plane  Containing  Apertun* 


Figure  47.  Descriptions  of  Aperture,  Magnetic  Field  and  Point  Under 
Considerations:  a)  Looking  Down  on  XZ  Plane  b)  Looking 
From  Outside  Onto  XY  Plane 


aperture  is  assumed  to  lie  parallel  to  the  plane  defined  by  the  reference  X and 
Y axes.  Frequently  the  XYZ  zero  point  of  the- reference  axis  will  be  taken  to 
coincide  with  the  center  of  the  aperture.  While  this  is  convenient,  it  if  not  nec- 
essary. The  aperture  is  assumed  to  be  an  ellipse  with  a major  axis,  LI,  and  a 
minor  axis,  L2.  The  angle  which  the  major  axis  makes  with  respect  to  the 
reference  axis  is  called  ANGA. 

The  coordinates  of  the  center  of  the  aperture  — XA,  YA,  and  ZA  — are 
input  quantities.  Likewise,  the  length  of  the  major  axis,  LI.  and  the  length  of 
the  minor  axis,  L2,  and  ANGA  are  also  input  quantities. 

The  magnetic  field  which  illuminates  the  aperture  is  assumed  to  lie  in  the 
same  XY  plane  as  that  containing  the  aperture.  The  field  vector  is  oriented 
at  an  angle  ANGH  to  the  reference  X axis.  The  magnitude  of  the  external  mag- 
netic field,  HEXT,  and  the  angle  it  makes  with  respect  to  the  reference  X axis, 
ANGH,  are  also  input  quantities. 

A defined  quantity  used  during  the  running  of  the  program  is  the  angle  be- 
tween the  major  axis  of  the  aperture  and  the  magnetic  field  vector,  ANAH. 

The  point  under  consideration  is  defined  in  terms  of  the  reference  X,  Y, 
and  Z axes  by  the  parameters  XP1,  YP1,  and  ZP1.  These  are  not  input  quan- 
tities. During  the  running  of  the  program,  the  coordinates  defining  the  pomt 
are  translated  to  a new  set  of  axes,  defined  by  the  major  and  minor  axes  of 
the  elliptical  aperture.  These  latter  are  not  shown  on  Figure  47,  but  go  by 
the  designations  XP2,  YP2,  and  ZP2. 

A detailed  flow  chart  of  the  MAIN  program  is  given  on  Figure  48.  The 
program  starts  with  a series  of  comments  on  the  program  and  a set  of  abbre- 
viated operating  instructions.  These  are  given  in  lines  100-1600.  They,  of 
course,  do  not  affect  the  running  of  the  program.  Definitions  and  dimensions 
of  file  names,  variables,  and  arrays,  are  given  in  lines  1610-1650. 

The  name  of  the  file  holding  the  input  data  is  given  in  lines  1660  and  1680. 
Since  the  program  is  at  present  configured  for  the  Geneial  Electric  time 
sharing  system,  this  name  is  an  input  quantity  entered  on  the  teletypewriter. 

For  batch  processing  a change  will  have  to  be  made  at  this  point. 

Lines  1700-1960  are  devoted  to  housekeeping  and  the  setting  up  of  for- 
mats of  headings. 

The  coordinates  defining  the  center  of  the  aperture  are  read  at  line  1970.  The 
date  on  the  size  and  orientation  of  the  aperture  is  read  at  line  1 980,  and  data  on  the  mag- 
nitude and  orientation  of  the  external  magnetic  field  is  read  at  line  2050. 

The  next  quantity  read  is  Dl,  a flag  used  to  say  whether  or  not  there  is  a 
reflecting  surface  to  be  considered.  D2  defines  the  location  of  this  reflec- 
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Figure  48.  MAIN  Program  — Detailed  Flowchart  (Sheet  1 of  5) 
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Figure  48.  MAIN  Program  --  Detailed  Flowchart  (Sheet  2 of  5) 
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Figure  48.  MAIN  Program  — Detailed  Flowchart  (Sheet  3 of  5) 


Figure  48. 


i i nes  33'H}-3620 


1 inc  .'iG.'iO 


Line  3650 


l ine  3670 


Lim  .1600 


DO  statements 
at  lines  3380 
and  3300 


Lines  3750-3770 


DO  statement 
at  line  3710 


Line  3780 


Lines  3800-3880 


MAIN  Program  — Detailed  Flowchart  (Sheet  4 of  5) 
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ting  surface  along  the  axis.  If  a reflecting  surface  is  not  to  be  considered,  a 
dummy  value  is  read  at  this  point.  Depending  on  the  value  of  D1  the  appro- 
priate heading  is  then  printed.  The  next  quantities  read  are  D3  and  D4.  D3  is 
the  flag  used  to  indicate  whether  or  not  a tabulation  is  desired  of  the  fields 
behind  the  aperture.  D4  describes  whether  the  tabulation  of  field  intensities 
are  to  be  printed  in  rectangular  or  spherical  coordinates.  A dummy  value 
must  be  entered  even  if  the  fields  are  not  to  be  tabulated.  Next  read  are  nine 
quantities  defining  the  volume  over  which  the  magnetic  fields  are  to  be  tabulated. 
These  are  used  to  set  up  the  appropriate  DO  loops. 

ZPA  defines  the  point  along  the  Z axis  at  which  the  tabulation  is  to  start, 
and  ZPB  defines  the  point  at  which  the  tabulation  is  to  stop.  ZPC  defines  the 
interval.  YPA,  YPB ZBC  define  similar  quantities  of  the  X and  Y axis. 

The  effective  dipole  moments  depend  upon  the  size  of  the  aperture.  The 
quantities  which  control  these  dipole  moments,  All  and  A22,  are  calculated 
using  the  subroutine  SHAPEFAC.  This  subroutine  is  a straightforward  evalu- 
ation of  the  equations  given  earlier  in  this  section  under  "Theory,  " and  so  is 
not  further  described  by  flow  charts.  The  quantity,  A33  which  is  also  evaluated 
by  SHAPEFAC,  is  not  used  in  this  program.  It  relates  to  the  effective  elec- 
tric dipole  moment  if  there  were  an  electric  field  impinging  on  the  aperture. 
While  the  computation  routines  and  housekeeping  routines  to  be  described  later 
would  evaluate  the  effects  of  an  electric  field,  they  have  not  been  incorporated 
in  this  program  at  this  time. 

Lines  2250  through  2400  relate  to  housekeeping  and  are  self-explanatory 
in  Figure  48.  At  the  end  of  this  housekeeping,  there  will  have  been 
generated  a set  of  coordinates  of  the  point  at  which  the  magnetic  field  is  de- 
sired. This  magnetic  field  is  calculated  with  the  subroutine  MAGFLD,  which 
is  described  below.  The  magnetic  field  strengths  returned  by  MAGFLD  are 
then  printed  in  either  rectangular  or  spherical  coordinates  as  requested  by 
the  input  data.  The  process  by  which  rectangular  coordinates  are  resolved 
into  spherical  coordinates  is  given  on  a separate  detailed  flow  chart. 

When  the  above  calculations  have  been  finished,  the  quantity  D5  is  read. 
This  quantity  is  a flag  used  to  indicate  whether  or  not  calculations  are  required 
of  the  flux  through  a defined  loop.  If  these  calculations  are  not  required,  the 
program  stops.  If  they  are  required,  a counter,  JX,  is  set  to  zero  and  the 
XYZ  coordinates  of  the  4 points  outlining  the  desired  loop  are  read.  These 
steps  occupy  lines  2860  through  2890. 

In  lines  2920  through  3170  are  calculated  the  lengths  of  the  sides  and  di- 
agonals of  the  quadrilateral  loop,  and  from  them  the  area  of  the  loop. 

In  lines  3190  to  3280,  the  distance  horizontally  across  the  loop  at  its 
midpoint  is  calculated  "Horizontal"  is  here  taken  to  be  the  direction  from 
point  P2  to  Point  P3  or  from  point  PI  to  point  P4 . The  term  "vertically"  is  taken 


114 


to  be  in  the  directionfrom  point  PI  to  point  P2  or  point  P4  to  point  P3 . These  terms 
in  this  sense,  have  no  relation  to  whether  the  loop  itself  is  oriented  horizon- 
tally or  vertically  with  respect  to  the  reference  XYZ  axes. 

In  lines  3300  through  3370  are  calculated  the  components  of  the  unit  vec- 
tor perpendicular  to  the  plane  defined  by  the  loop  under  consideration.  Mathe- 
matically this  operation  consists  of  taking  two  vectors  that  lie  in  the  planes 
point  1 - point  2,  and  point  2 - point  3,  and  taking  the  cross-product  of  these 
two  vectors. 

Next,  the  quadrilateral  is  divided  into  twelve  strips  vertically  and  twelve 
strips  horizontally,  and  the  field  strength  calculated  at  the  intersection  of 
each  of  the  dividing  lines.  This  makes  a total  of  169  points.  This  field 
strength  is  calculated  by  a call  to  the  subroutine  MAGFLD.  MAGFLD  returns 
the  X,  Y,  Z components  of  field  strength  with  reference  to  the  original  refer- 
ence axes. 

In  line  3650,  a dot  product  of  the  field  strength  vector  and  the  unit  vector 
normal  to  the  plane  is  performed  in  order  to  determine  the  component  of  the 
magnetic  field  perpendicular  to  the  loop  under  consideration.  These  field 
strengths  are  loaded  into  an  array,  HN,  at  line  3670.  The  vertical  distances 
along  each  of  the  thirteen  strips  are  loaded  into  an  array  T86  at  line  3690. 

In  lines  3750  through  3770,  line  integrals  of  the  magnetic  field  strength 
are  taken  along  each  of  the  thirteen  vertical  paths.  This  integral  is  evaluated 
numerically  by  dividing  the  vertica1  strip  into  two  sections  of  seven  points  each. 
One  point  is  common  to  each  of  the  two  sections.  The  integration  routine  used 
is  Weddle's  rule,  which  was  described  under  "Theory. " These  line  integrals 
are  stored  in  the  array  PATHA.  The  thirteen  line  integrals  are  then  integrated 
horizontally  to  obtain  the  total  magnetic  field  linking  the  plane.  The  double  in- 
tegral of  H is  taken  in  line  3860  and  then  multiplied  by  the  permeability  of  air 
to  obtain  the  total  magnetic  flux  in  webe.vs.  This  latter  multiplication  is  taken 
at  line  3870. 

Finally,  the  loop  number,  JX,  the  coordinates  of  the  points  defining  the 
loop,  the  area  of  the  loop,  and  the  total  magnetic  flux  linking  the  loop  are 
printed  out.  The  program  then  loops  back  to  read  in  the  coordinates  of  ad- 
ditional points,  if  there  are  any  additional  loops  to  be  considered.  If  no  data 
are  found  in  the  input  files,  the  program  stops. 

Rectangular  to  Spherical  Coordinates 

Figure  49  is  a detailed  flow  chart  of  the  process  by  which  the  rectangular 
coordinates  of  the  magnetic  field  are  resolved  into  spherical  coordinates. 

Figure  50  shows  the  conventions  regarding  the  designation  of  the  spherical 
coordinates.  There  are  two  angles  to  be  calculated:  the  latitude  angle,  the 
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HPXl  * Field  strength  in  X Direction 
'.il’Yl  = Field  strength  in  V Direction 
HPZ1  « Field  atrength  In  Z Direction 


Ang  1 = ^^jpYl) 


Ang  1 * 90-tar 


Angle  1 may  be  in  either  Qiadrant  I or 
Quadrant  IV. 


PPX1MD| 


Ang  2 * tan"1 


(4030)  iL  Angle  2 may  be  in  any  quadrant. 

^ HPY1  \_!i2 

V positive^  1 


(4040)  y 
(4120)  1 Yea 


1 (4050) 


Angle  1 ia  in  Quadrant  IV 
Ang  1 - 180- ang  1 


(4120)  Line  2680 


Angle  may  be  in 
either  Quadrant  1 
or  Quadrant  IV. 


la 

HPZ1 

poaitive? 


la 

HPZ 

poaitive? 


Angle  2 ia  in  Quadrant  IV 
ang  2 - -ang  2 


» Yea  (4190)  l ine  2780 


Angie  2 ia  in  Quadrant  1 
ang  2 11  ang  2 


iLine  2770 
(4210) 


Angle  2 ia  Quadrant  111 
ang  2 - -(160- ang  2) 


4215  Line  2780 


Figure  49.  Detailed  Flowchart 
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Figure  50.  Conventions  Regarding  Angles:  a)  Latitude  and  Longitude 
Angles  Defined;  b)  Quadrant  Designations  for  Longitude; 
c)  Quadrant  Designations  for  Latitude 


angle  with  respeet  to  the  vertical  Y axis;  and  the  longitude  angle,  the  angle 
in  the  XZ  plane  from  the  positive  X axis.  Lest  there  be  confusion  as  to  why 
the  positive  X axis  points  to  the  left,  remember  that  Figure  50  shows  the  re- 
gion behind  the  aperture.  When  viewed  from  outside,  where  the  magnetic 
field  originates,  the  reference  X axis  has  its  positive  sense  to  the  right.  The 
latitude  and  longitude  angles  are  e ailed  ANG1  and  ANG2,  respectively,  in  the 
program. 

These  angles  are  basieally  calculated  from  the  arc  tangent  of  the  respec- 
tive components,  HX  and  HZ  for  ANG2  (LONG)  and  HY  and  HD  for  ANG1  (LAT). 

D is  the  length  of  the  projection  of  the  H field  veetor  in  the  XZ  plane.  There 
are  two  problems  in  this  resolution.  The  first  is  to  ensure  that  under  no  condi- 
tion does  the  denominator  in  the  argument  for  the  are  tangent  go  to  zero.  If 
it  does  go  to  zero,  appropriate  angles  are  calculated,  but  annoying  error  mes- 
sages are  still  generated  and  printed  by  the  computer.  This  is  prevented  from 
occurring  by  the  switeh  at  statement  4020,  line  2630,  and  the  alternate  methods 
of  calculating  the  angle  at  statements  4010  and  4004.  The  appropriate  switch 
and  statements  for  ANG2  oeeur  at  statements  4012,  4014,  and  4020. 

The  second  problem  relates  to  determining  the  appropriate  quadrant  in 
which  the  angle  Ues,  since  the  arc  tangent  routine  does  not  intrinsically  re- 
solve quadrants.  Upon  evaluation  of  the  arc  tangents,  angle  1 may  be  in  either 
quadrant  1 or  4.  Quadxants  2 and  3 are  forbidden  regions,  beeause  the  polarity 
of  the  D component  is  always  positive,  inasmuch  as  it  is  taken  by  the  veetor 
addition  of  the  X and  Z components.  Angle  2 may  be  in  any  quadrant.  The 
switches  at  statements  4030,  4120,  4140,  and  4.160,  resolve  the  question  of 
appropriate  quadrants.  Appropriate  statements  add  or  subtract  180°  or  re- 
verse the  sign  of  the  angles.  The  logic  is  straightforward,  though  a bit  in- 
volved, and  is  shown  on  the  remainder  of  Figure  49. 

MAGFLD  Program 

The  major  subroutine  used  in  the  program  APERATURE  is  MAGFLD. 

Figure  51  is  the  flow  chart  for  this  subroutine.  The  program  is  entered  at 
line  4320,  using  the  quantities  shown  at  the  top  of  Sheet  1.  After  the  initial 
quantities  are  defined,  the  first  task  performed  is  to  translate  the  coordinates 
of  the  point  under  study  from  the  original  reference  X and  Y axes  to  a new  set 
of  axes,  oriented  along  the  major  and  minor  axes  of  the  aperture.  This  is  done 
at  lines  4380  through  4430.  The  Z coordinate  of  the  point  under  study  is  also 
shifted  to  a new  Z axis,  centered  on  the  aperture.  The  distances  from  the  middle 
of  the  aperture  to  the  point  under  study,  in  terms  of  the  new  coordinate  geometry, 
are  then  calculated  in  lines  4480  through  4510. 

The  quantities  Cl  and  C2  are  then  calculated.  Cl  and  C2  are  basieally 
the  distances  from  the  center  of  the  aperture  to  the  point  under  study,  although 
they  also  inelude  a term  related  to  the  length  of  the  major  and  minor  axes  of 
the  aperture.  Accordingly,  these  terms  eannot  go  to  zero,  even  if  the  point 
under  study  were  to  be  at  the  center  of  the  aperture. 
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Quantities  Used  to  Enter  MAGFLD  Are: 


y 


ANGA  - Angle  aperture  makes  to  X axis 
ANGH  - Angle  field  makes  to  X axis 


} 

} 


XP1 
YP1 
ZP1 
XA 
YA 
ZA 
HEXT  - 
All  - 
A22  - 

LI 
L2 
D1 
D2 


X,  V,  Z coordinates  of  point  at 
which  field  is  to  be  calculated. 

X,  Y,  Z coordinates  of  center 
of  aperture 

External  magnetic  field  strength 
Shapefactor  for  major  axis  of  aperture 
Shapefactor  for  minor  axis  of  aperture 
Length  of  major  axis  of  aperture 
Length  of  minor  axis  of  aperture 
YES/NO  as  regards  reflecting  surface 
Z coordinate  of  reflecting  surface 


Quantities  Returned  by  MAGFLD  Are: 


HPXn 
HPY 1 i 
HPZ1  J 


Magnetic  field  strengths  in  X,  Y.  Z 
directions  at  point  under  calculation 


Lines  4340-4370 


Lines  4380-4430 


Figure  51.  Flow  Chart  for  Subroutine  MAGFj_,D  (Sheet  1 of  3) 
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Figure  51.  Flow  Chart  for  Subroutine  MAGFLD  (Sheet  2 of  3) 
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(9140) 


Calculate  F 8. 


Lines  4790*4800 


(9970) 

/VeryX 

,Xclose  to  the 
aperture  In  this  Y^1 
X direction’/"^ 


Unes  4830*4830 


Calculate  P 2. 


Lines  4850*4880 


Calculate  the  field  strength  with 
respect  to  shifted  coordinates. 
(HPX).  (HPV),  QIP2) 


Lines  4910*5010 


'•  \ 

* a reflecting 

surface?  y/ 


Yes  (D1  - 1) 
0682) 


Add  new  field  strengths  due  to  ap- 
erture  to  those  stored  in  registers 
HPX2 
HPY2 


Lines  5030  5050 


Calculate  the  position 
of  the  next  aperture 
Image. 


X Was  the  X 
y'  imount  added  slgnlf-^^ 
leant  compared  to  what  was 

\ already  In  the 
registers  ? 


Lines  5060*5080 
Statement b 9690 
No  ®700 
— | 9710 


i - n*  i i ■ 

U Lint-  111 


X Have  >. 
ten  images 

been  considered  yet? 
N is  the  S 
^Xcounter.V^ 


Add  new  field  strengths  due  to  ap- 
erture to  tho*;e  stored  in  registers 
HPX2 

] HPY2 

11PZ2 


Unes  5130-5150 


Cslcuiste  field  strengths  with 
reference  to  the  origtnsl  XYZ 
axis. 


Lin.  3160-5200 


Figure  51.  Flow  Chart  for  Subroutine  MAGFLD  (Sheet  3 of  3) 
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In  subsection  "Theory,  " the  field  strength  at  the  point  under  study  is 
seemed  to  be  that  due  to  the  magnetic  moment  of  the  dipoles  formed  by  the 
major  and  the  minor  axes  of  the  elliptical  aperture.  These  dipole  moments 
are  the  product  of  the  magnetic  field  strength  along  the  major  and  minor  axes 
of  the  aperture,  the  lengths  of  the  major  and  minor  axes,  and  the  shape  fac- 
tors for  the  aperture.  These  shape  factors,  which  also  include  the  lengths  of 
the  major  and  minor  axes,  were  calculated  in  the  subroutine  SHAPEFAC. 
These  factors  are  all  calculated  in  lines  4580  through  4640. 

In  subsection  "Theory,  " die  field  strength  ev.  ations  were  presented  in 
terms  of  the  subfactors  FI  through  F4  and  G1  through  G8.  In  die  subroutine, 
the  quantities  FI  and  F2  are  the  same  as  the  quantities  FI  and  F2  derived  in 
the  subsection  on  theory;  quantities  F3  and  F4  mentioned  in  there  are,  how- 
ever, replaced  by  their  corresponding  equivalents,  F7  and  F8.  The  quanti- 
ties F3  through  F6  in  the  subroutine  are  not  related  to  any  corresponding 
quantities  in  the  subsection  on  theory. 

The  quantities  FI,  F2,  F7  and  F8  involve  raising  the  quantities  Cl  and 
C2  to  powers  up  to  and  including  the  7th  power.  When  the  point  under  study 
is  very  close  to  the  X and  Y axes,  but  not  exactly  on  the  axes,  underflow  con- 
ditions are  generated  in  the  computer.  Correct  numerical  answers  are  re- 
turned, but  annoying  error  messages  are  still  printed.  In  order  to  eliminate 
these  error  messages,  there  are  switches  at  statements  9165,  9195,  9230, 
and  9270  which,  when  appropriate,  calculate  the  quantities  FI,  F2,  F7,  and 
F8  by  their  small  argument  equivalents.  This  process  of  evaluating  the  terms 
in  the  field  equations  occupies  the  space  from  lines  4650  through  4860.  The 
components  of  the  magnetic  field  strengths  are  then  evaluated  at  lines  4990 
through  5010. 

At  this  point  the  presence  or  absence  of  a reflecting  surface  must  be 
treated.  If  a reflecting  surface  is  not  present,  as  indicated  by  the  switch  at 
line  5020  or  in  statement  9680,  the  calculated  magnetic  field  vectors  are  ro- 
tated back  to  the  original  reference  axis  in  lines  5180  through  5200  and  the 
quantities  HPX1,  HPY1,  and  HPZ1  are  returned  to  the  program  MAIN. 

If  a reflecting  surface  is  present,  the  field  components  calculated  are 
added  to  the  contents  of  the  storage  registers  HPX2,  HPY2,  and  HPZ2.  (Ini- 
tially these  storage  registers  had  been  set  to  zero  at  lines  4440  through  4460. ) 
The  program  then  determines  whether  the  quantities  added  to  the  storage 
registers  HPX2,  HPY2,  and  HPZ2  were  significant  compared  to  what  was  al- 
ready in  the  registers.  For  this  first  loop  through  the  program  the  quantities 
of  course  were  significant;  there  was  nothing  stored  in  those  registers  to  be- 
gin with.  The  program  then  loops  back  to  calculate  the  field  strengths  pro- 
duced by  the  first  reflection  of  the  aperture  in  the  reflecting  surface.  The 
position  of  the  reflection  is  along  the  Z axis  at  a spacing  from  the  original 
aperture  equal  to  twice  the  spacing  to  the  reflecting  surface.  This  new  posi- 
tion along  the  Z axis  is  calculated  at  line  5110,  statement  9740. 
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The  program  then  adds  the  field  strengths  produced  by  successive  reflec- 
tions to  those  stored  in  the  register*?  HPX2,  HPY2,  and  HPZ2,  testing  at  each 
time  to  see  whether  the  contribution  from  the  aperture  under  study  was  signi- 
ficant enough  to  bother  with.  The  number  of  times  through  this  loop  is  counted 
with  the  counter  N.  The  coordinate  of  the  aperture  under  study  increases 
rapidly  as  the  program  goes  through  this  cycle,  and  eventually  the  contribution 
to  the  total  field  strength  from  the  higher  order  reflections  becomes  negligible. 
The  counter  N is  used  at  the  switch  point  9720  to  break  out  of  the  loop  if  the 
field  strength  has  not  converged  to  its  final  value  after  treating  ten  images. 

If  the  contribution  from  the  last  image  was  negligible,  or  if  ten  images  had 
been  considered,  the  field  strengths  in  the  registers  are  then  rotated  back  to 
the  original  X,  Y,  Z axes,  and  the  quantities  HPX1,  HPY1,  and  HPZ1  are  re- 
turned to  the  program  at  MAIN. 

To  ease  the  task  of  going  through  the  program.  Table  1 lists  the  statement 
numbers  in  ascending  order  versus  their  Corresponding  line  numbers.  If  the 
program  were  to  be  sequenced  the  line  numbers  would  change. 

VAUDATKDN  OF  APERTURE 


Validation  of  the  computer  program  APERTURE  is  based  upon  a compari 
son  of  the  computer-generated  results  with  the  results  predicted  by  classical 
electromagnetic  theory  — that  the  flux  density  due  to  a magnetic  dipole 
decreases  as  a function  of  l/r3  for  large  values  of  r.  The  computer  results 
are  shown  in  Figures  52  through  54. 

These  figures  show  the  computer  results  to  be  in  agreement  with  this 
l/r3  decrease.  The  orientation  of  the  vector  field  is  also  equivalent  to  that 
predicted. 
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Table  1 


STATEMENT  NUMBERS  VERSUS  LINE  NUMBERS 


Line  No. 

Statement  No. 

Line  No. 

Statement  No. 

10 

1660 

4000 

2550 

IS 

1670 

4002 

2560 

20 

1660 

4004 

2570 

SO 

1690 

4006 

2580 

110 

1710 

4010 

2590 

115 

1720 

4012 

2600 

120 

1730 

4014 

2610 

122 

1740 

4016 

2620 

123 

1750 

4020 

2830 

130 

1770 

4030 

2640 

140 

1790 

4040 

2650 

145 

1800 

4050 

2660 

ISO 

1610 

4110 

2670 

155 

1820 

4120 

2660 

160 

1630 

4130 

2690 

165 

1640 

4140 

2700 

170 

1850 

4150 

27  0 

175 

1860 

4160 

2720 

180 

1670 

4170 

2730 

185 

1860 

4180 

2740 

186 

1900 

4190 

2750 

180 

1910 

4200 

2760 

192 

1920 

4210 

2770 

195 

1930 

4215 

2760 

200 

1940 

4230 

2600 

210 

1950 

4315 

3890 

220 

1960 

4320 

4000 

275 

2110 

4330 

4010 

260 

2120 

4340 

4020 

265 

2140 

4380 

4040 

290 

2150 

5230 

5220 

295 

2170 

9130 

4490 

1002 

2470 

9140 

4500 

1199 

2460 

9150 

4510 

1210 

2520 

9160 

4520 

1400 

9165 

4860 

1750 

2490 

9170 

4670 

1950 

2610 

9160 

4700 

1955 

2150 

9190 

4710 

1957 

2670 

9195 

4720 

1960 

4050 

9200 

4730 

2840 

9210 

4760 

2280 

9220 

4770 

2201 

2270 

9230 

4780 

2202 

2280 

9240 

4790 

2204 

2360 

9250 

4020 

2206 

2310 

9260 

4830 

2207 

2320 

9270 

4640 

2206 

2340 

9260 

4650 

2450 

2370 

9300 

4680 

3130 

4280 

9320 

4900 

3140 

4270 

9660 

5020 

3160 

4290 

9862 

5030 

3170 

4300 

9890 

5060 

3370 

2510 

9700 

5070 

3550 

3560 

9710 

5060 

3560 

3390 

9720 

5090 

3570 

3400 

9730 

5100 

3860 

3700 

9740 

5110 

3690 

3710 

97  CO 

5130 

3950 

3790 

9770 

5155 
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Figure  53.  Field  Intensity  - Y Component 


a 


126 


Section  4 


CONCLUSIONS 


Two  computer  programs,  APERTURE  and  DIFFUSION,  have  been  developed 
for  the  calculation  of  probable  electromagnetic  fields  and  resulting  induced 
voltages  in  aircraft  electrical  circuits.  These  programs  should  enable  the  re- 
searcher or  system  designer  to  determine  the  order  of  magnitude  of  lightning 
induced  voltages  which  may  be  induced  in  simple  circuit  geometries  by  lightning 
strokes  of  any  assumed  amplitude  and  waveform.  By  variation  of  the  conductor 
location  parameters,  the  programs  enable  the  designer  to  determine  the  best 
location  (i.e. , where  coupling  is  minimized)  within  the  airframe  for  placement 
of  conductors. 

The  DIFFUSION  program  is  based  on  calculation  of  the  magnetic  fields 
which  occur  inside  an  airframe  as  a result  of  lightning  current  diffusing  to 
the  inside  surface  of  its  metallic  skin.  The  program  therefore  assumes  that 
the  airframe  skin  is  metallic  and  has  no  apertures.  This  is  the  flux  which 
normally  exists  inside  an  all-metallic  airframe,  and  should  be  considered  as 
the  minimum  to  which  internal  fields  can  be  reduced  in  a metallic  airframe  of 
given  skin  material  and  thickness  by  such  means  as  closure  of  apertures  and 
improvements  in  electrical  bonding. 

Because  diffusion  fields  are  of  relatively  low  amplitude  and  slower  rates 
of  rise  than  their  external  counterparts,  voltages  induced  by  diffusion  fields 
linking  small  circuit  loop  areas  such  as  those  formed  between  parallel  pair 
or  twisted  pair  conductors  are  likely  to  be  small.  On  the  other  hand,  large 
loops,  such  as  those  formed  between  either  conductor  of  a pair  and  the  airframe, 
may  receive  high  induced  voltages  from  diffusion  fields.  This  is  especially  true 
because  the  diffusion  fields  are  usually  present  throughout  the  entire  length  of 
such  a circuit. 

The  APERTURE  program  calculates  the  fields  penetrating  the  interior  of 
the  airframe  from  a given  field  tangential  (in  any  assumed  direction)  to  the 
outside  surface  of  the  airframe  at  the  aperture  in  question.  These  fields  pene- 
trate directly  into  the  interior  of  the  airframe  but  are  strong  only  in  the  vicinity 
of  the  aperture.  If  a parallel  pair  of  conductors  passes  nearby,  the  aperture 
fields  are  often  of  great  enough  amplitude  and  rate  of  rise  to  induce  large  volt- 
ages. If  located  some  distance  away  from  the  aperture,  however,  resulting  in- 
duced voltages  may  be  small,  because  the  field  intensity  falls  off  as  the  square 
or  cube  of  the  distance  from  the  aperture. 


Thus  a complete  analysis  of  a particular  situation  will  usually  require  the 
use  cf  both  computer  programs  and  superposition  of  the  results  of  one  on  those 
of  the  other  for  consideration  of  the  worse  case. 
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At  present,  APERTURE  and  DIFFUSION  deal  with  relatively  basic  geom- 
etries and  do  not  account  for  such  details  as  internal  structural  components 
(e.g. , spars  and  ribs;,  concentration  of  lightning  current  around  the  points 
of  stroke  entry,  or  leakage  through  resistive  joints  or  bonds.  It  will  therefore 
be  desirable  to  develop  further  refinements  to  permit  consideration  and  ac- 
curate calculation  of  the  effects  of  such  details  as  ribs,  spars,  seams,  access 
doors,  flap  openings,  as  well  as  such  other  objects  as  antennas  and  radomes. 
Each  of  these  additions  should  be  validated  by  comparison  with  measured  test 
data  obtained  from  other  programs  of  aircraft  lightning  induced  voltage  mea- 
surement. 

It  may  also  be  advantageous  to  convert  the  input  and  output  formats  of  the 
programs  to  the  same  format  as  the  one  used  in  the  Air  Force  intersystern 
analysis  program  (IAP).  The  latter  is  a frequency -domain  input/ output  format 
which  expresses  intersystem  electromagnetic  interference  (EMI)  in  terms  of 
its  frequency  spectral  content  (energy  at  each  frequency  within  a wide  bandwidth 
of  frequencies).  Basically,  conversion  of  the  basic  lightning  induced  voltage 
model  to  this  format  will  require  conversions  of  the  calculated  induced  impulse 
voltages  to  their  Fourier  spectral  coefficient  equivalents;  the  interference  from 
lightning  is  therefore  expressed  in  the  same  frequency  spectral  language  as  the 
EMI  already  calculated  by  the  IAP.  Other  format  changes  will  involve  airframe 
geometrical  descriptions.  These  changes  are  not  expected  to  be  extensive, 
however,  and  if  made  may  promote  use  of  these  lightning  induced  voltage  models 
by  engineers  concerned  with  the  solution  of  related  EMI  problems  as  well. 
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Appendix  i 

DERIVATION  OF  TOTAL  FLUX 


The  purpose  of  the  derivation  given  here  is  to  find  the  total  flux,  f,  gen- 
erated by  the  current  from  a given  filament,  passing  through  an  area  bounded 
by  lines  parallel  to  the  wire  at  distances  D1  and  Da  from  it,  along  the  wire 
from  A 8 to  i r Points  A t and  A 3 (shown  in  Figure  55)  are  the  beginning  and 
end  of  a circuit  conductor. 


Figure  55 . Geometry  of  Flux  Derivation 


This  comprises  three  separate  integrals: 
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.*£  fL' 

4*  D, 


a + r* 


dr 


4- 


fDS  r ^ 

— d°g  (21  - 2L  + 2 Ve  8 - 2 Lj f + La  + r3) ) 
4*  JD  Lr 


^ 2 


dr 


^ l 


,Da  f 

if°L  l-L  ( Vi  8 “ 2Li  4-  L*+  ra  + L log  (2j 8 - 2L 

4,r  JD,  lr 


+ 2 V*  a - 2 Li  + La+  r8 
Rearranging  terms  and  substituting  i l and  j ? s> 

- V*ia  + rg 


dr 


•])  - ^ 
47T 


ri* 

JD, 


2a+  r8 


+ Vl,a  - 2 Li  , + L8  4-  re,  - Via8  - 2Ii,4  L8  4- 


(125) 


4-  ~ log  (21  a - 2L  4-  2 Vis3  - 2Lia4  Ls+  r8) 


— log  (2ia  - 2L  4-  2 Vi88  - 2Li2  4-  L°+  r8 


4-  ~ log  (2ij  - 2L  + 2 Vii8  - 2Lix  + L8+  r8) 


- ~“log  <2lx  - 2L  4-  2 Vi8  - 2Lix  4-  L8+  r2) 


dr 


(126) 
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dr 
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Appendix  It 
FUSELAGE 


The  geometry  of  the  fuselage  model  is  broken  down  to  straight  line  segments 
and  circular  sections.  These  configurations  can  be  used  to  describe  the  front 
view  of  most  fuselage  geometries.  To  use  the  program  the  fuselage  is  laid 
out  as  in  Figure  56.  The  left  side  and  bottom  are  along  the  XY  axis.  The  top 
is  at  Y =Yj  and  the  right  side  at  X =X1.  The  circular  sections  in  the  corners  are 
set  up  so  that  the  radii  of  the  top  sections  are  the  same  and  the  lower  sections 
are  the  same. 


Figure  56.  Fuselage  Front  View 

It  is  anticipated  that  typical  aircraft  will  be  symmetrical  from  side  to  side 
on  the  top  and  on  the  bottom,  but  that  the  lower  section  may  differ  from  the 
upper  section  in  curvature.  This  method  can  be  used  directly  on  a structure 
that  has  no  straight  segments,  representing  the  minimum  program  to  encompass 
the  maximum  geometry  anticipated. 

If  the  fuselage  is  tapered,  the  program,  which  assumes  no  taper,  can  be 
sectioned  into  two  or  more  straight  pieces;  the  same  setup  would  be  used,  with 
different  physical  sizes. 

The  front  view  of  a typical  fuselage  is  shown  in  Figure  56.  The  distances 
Dj  and  Ds  are  calculated  for  the  filaments  along  the  straight  line  sections  and 
the  four  curved  sections. 

In  the  straight  sections  of  Figure  56, 

Dj  = 'JoT  X.)»V<T-  Y,)B  ’ (129) 

• For  the  section  from  0,  k3  to  0,  k2: 

X = 0,  step  Y in  J steps  from  k.  to  k2 
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• For  the  section  from  hg,  to  hjY^ 
Y=Ylf  step  X in  J steps  from  kg  to  h, 

• For  the  section  from  X^kj  to  Xjk4:- 

X = Xj,  step  Y in  J steps  from  kj  to  k4 

• For  the  section  from  h4,  0 to  ha  0: 

Y=0,  step  X in  J steps  from  ha  to  h4 

For  the  curved  sections  use  Figure  57 . 


Figure  57  . Curved  Sections  Geometry 

For  the  first  and  second  quadrants, 

kj  = k8  = k 
sin  0 - y."k 

Ej 

y = Rj  sin  6 + k 
The  equation  for  each  section  is; 

(X  -h) 8 + (Y-k)8  = Rj8  (130) 

with  h and  k as  h^  for  the  first  quadrant  and  kaha  for  the  second  quadrant: 

(X  - h)8  = R 8 - (Y  - k)8 

X = VR*  - iY-k)8  + h 
Substitute  Y = Ra  sin  6 + k: 

X = (Rj  sin  e + k - k)8  + h = VRf  " (Rj  sin  6)a  + h 

Di  = V(X  - X,)3+  (y  -yA)s 


(131) 

(132) 
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(133) 

(134) 


Substitute  for  X and  Y 


Di  = '\/VR1a  - (Ri  sin  0 )a  + h - XA)8  + (R1  sin  6 + k - YA)S 


(135) 


Da  is  Dx  with  XbY8  substituted  for  XAYA: 

Dx  = ^( Rj  cos  0 + h - XA)a+  (R1  sin  9 = k - YA)S 

First  quadrant:  step  9 from  0 to  ^ in  — steps. 

i K^ 

Second  quadrant:  step  9 from to  ir  in  — increments. 

4 K, 


(136) 


For  the  third  and  fourth  quadrants: 

k-Y 


sin  9 = 


R 


Y = k-R  sin  9 


X is  still  the  same,  with  Ra  substituted  for  Rji 


X = VR®  - (Ra  sin©)3  -t-h 


D1  = V(X-Xa)8+  (Y  -Y  A)a 


(137) 

(138) 


=~\j  \/(Raa  - (Ra  sin  0)8  + h - XA)Z+  (k  - RJ  sin  9 - Y A)2  (129) 


Da  is  Dj  with  XbYb  substituted  for  XAYA: 

Di  = V(R,  cos  0 + h - XA)a  ~(k  - Ra  sinS- y7)“ 

3 J 

Third  quadrant:  step  9 from  n to  — 7r  in—  increments. 

^ Kg 

3 J 

Fourth  quadrant:  step  0 from  — n to  2 ir  in  — increments. 

2 Kg 


(140) 
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Appendix  III 

PROGRAM  LISTINGS  K>£  CDC6600 COMPUTER 


DIFFUSION  and  APERTURE  programs  were  supplied  by  General  Electric 
Corporate  Research  and  Development  to  Wright-Patterson  Air  Force  Base  for 
use  on  a CDC6600  computer.  The  DIFFUSION  program  listing  for  the  CDC  - 
6600  is  shown  in  Figure  58,  the  APERTURE  listing  in  Figure  59. 
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Figure  58.  DIFFUSION  Program  Listing  (Sheet  1 of  9) 
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iit^lHi 
1 1 M2  291 
ilMflN 
1 1 §021X1 


1 1 SI  2121 
IIIS2SM 
Ilil2l4§ 
IIII2IH 

OlHtHi 


AC  NINO  1 
AF NINO  2 
RCNI  NO  S 
ftCHTNO  4 
ftCNTNO  5 

:~57rS * - 

IFf ICtC0.?7IG0T0  IISI 

X1«*9 

V9«VS 

U«LS 

L2*L6 

Trio\.tr.inc3904k0 

G0T066I 

42f  l9>K5*t0SMIC‘04n 
IF  ID^.LT.ICIGOTOWO 
440  IF(07.LC  .mCOTOAAfi 
G0T04X9 

T4«Ts*m^nr-xiTi) 
!F(0K,LT,!CI  GOTO990 
411  XFIOA.LC.mAOTOMI 

ooiflsxs  ..  . 

Ill  Ll*LJ*t0JMIC~06tl 
TF(09,LT •lOCOTOtil 
123  (Ft07,LC .ICIG0T096Q 
COTOAvO 

_9*l  L?«L4*t04MIC-07ll 
T7  luiibl  •ItlCuiuiH 
GOTObll 

971  X4«X9*<0  J*  Wl 
GOTO460 

991  V9«V9* (06*04 ) 

GOTObll 

Ti|-ii«iTrm»wi 

C0T092I 

433  L2*L*»*  06*011 
641  N1UI-I1 
*t«*l-ftl 
N2«ftl 

nt2«Ti-ffi  ' ‘ 

M J«ft  j 

<!•« 

H4«Il-*2 

K4«R2 

XF  t A. CQ. All GOT 0760 


mnm 

IIH2HI 
1 1 If  TIM 
••Of 2411 
10002410 
II II 2421 
0I0024M 
1 1 112441 
IIII24II 

•111246ft 
00102471 
II 002491 
00012901 
• •••2911 


• •••2990 
00012941 
•I  II 2971 
011029ft! 

• III  2991 


769  CONTINUE 


ftCNINO  1 
REMIND  2 


IIM741I 
10102020 
•III  2411 
1910204ft 
IIII269I 
lllfTOI 

II 002010 

01014011 

lloilfot 

mimo 

•100  2721 
•I002T90 
Mil  2741 
TUI 2791 
11112761 

nnim 

IIII27II 

II1I270I 

IIII2IH 

irmni 

MM  2121 
1111200 
IIII2I4I 
11112091 
MII2I6I 
Tf  *1071 
1*112110 
CIII2I9I 
IIII29II 
IIII29X0 
Mil  2121 

Timm 

II II 2161 
IIII2990 
III! 2961 
•0112171 
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tii 


•orniai  ,n 

HITRH  TfT 


ffl  CONTINUE 


QaC/Nl 
”I*T  *“ 

UjMUl 

JMMH' 


yiUMIIMTOMI 
Ml  CONTINUE 

it  immwrmCTwm 


•IMltM 
9 Mm  ii 
IflMMtl 
iiMsm 
mm  mi 

IIMMM 

IINNM 

MMitn 

• •HUM 

muni 

IIMIltl 


ttr 


*»*I 
IW*91*l 


nr 


00«?iioutt>i*ii* 

WWR 

rot  #o*i 

iPiT.u.Ti-niuinu 

UN1TEI*)  N*T 

Ml  CINTtNUE 

111  TM  , 

HNU*9»J 

tray  .cr.Ti* » nwrowi 

UNITE <11  Mitt 

Ml  Ifjr.EQ.TlI&OTOlIM 
T«T1 


mm  mi 

•Hit  I 
MMliM 
I Will  71 
IIII11M 


•IMltM 

HMitn 

lllllttl 

•msm 

•IMltM 


001 I SI I01M*1 • 1UN 

I-I6UN-1  •*— 

immnCTTOwr 


JtlTEftl  xfr 

mrcwriNur  _ 

tIM  KUN«N1M 

wrnmsmi'H* 


nr 


K-KOUW-l 

TT.nW! 

XfaNl-Nf  USITlI 

If  Ut«M.tt>60T°il3l 
MITuc*  Klaft 


lit!  CONTINUE 

lifT  CfllfTIWdl- 


C NCN  NEXT 

wpcum 


nr 


001  tllKOUMl  t NUN 
Tl*K»0/Nt 


in 


Vt«Kt*lt»COSlTll 

irmrcr^ti  comm 

NNITCIU  Xt.ft 

‘"mrconTTRui 

111  I CONTINUE 

c tew  nm 


m 


tii 


LUMNlal 

— OOmVUHMil'Lt# 

(.■10UN-1 

' ~TI5t»D/Nl 

It*HlHtl*COI  I Til 

HMIMPSIHIMI 


irr 


nr 


hi 


IMXt.lE.HlItOTOltU 
— iwrrmw.T* 

Itll  CONTINUE 
~lftl  CtiTIMT 
C NCH  IE  SET 

WF=im 

OOUtlNOUNal.Ni# 

ir-umif-i 

T1«N*0/*1 

ftatft-YPsrNmt 

Tt>Kt«NlaCOSITll 

irmoninwwnn  ~ 

NNITEftl  Xt«  Yt 
131V' 'CONTINUE 
1JP I CONTINUE 

If  (XI.LE.lt) COTOlMI 

irm.ig.iiuoToim  * - 

If  (XI.GE.NlltOTOHM 
Trm.cr.HiuoToiwi 
lit 0 If  (YI.LE.ItttOTOlMI 
nil  IMVI.tE.XllCOTOUM 

comiii 

III!  If (XI.EQ.TIIGOTOltll 
Hill  - ~ 

TlaTI 

COTOlfM 

im  flair 

XT»N1*(VQIT( (Nil t (Nl*f IfYII  •■tl  1 1 
TlaTI 

60  TO  1PM 

■■  Tnr  irrmr  ;ur  cdtouii 

1ISI  If  (TI.CC.XllCOTOltM 
' - —flail 


•IMltM 
MMOTI 
iiuatii 
■IMltM 
•••HIM 
IIM111I 
•llllltl 
•IMVIM 
•iiiiiM 
••MUM 
•iiiiiM 
MM  Ilf  I 
IMIVMI 
0111  VIM 
••••$*•• 
•I II VI II 
IMIVMI 
MIIMM 

••••mi 

MM  MM 
IIMMM 
MMMPI 
MilMII 
MMMM 
MM  VIM 
MMVI1I 

miviti 

•Mlllll 

MMMM 

MMMM 

•IMVIM 

• Ml  Ilf  I 
MU  JVM 

MMMM 

MMM1I 

MM  lit! 

M 11 

•iMmo 

IIIIMII 

IIMMM 

MMMfl 

MMMlft 

MMMM 

urn 

• Ml  If  II 
MMIft9 
lllllfll 
IMIlfM 
MlllfM 
iiifstn 
IMIlffl 
lllllfll 
lllllfll 
•••HIM 

miovzi 
ii  ii  mi 
mi  mi 
tin  mi 
•Ml  Mil 
MUllfl 
•••HIM 
•IMVIM 
MM  II II 
Mimu 
IMIVMI 

•mini 

••MUM 

•Mlllll 

•IMVIM 

lllllttl 

•IIIIIM 

IMHIM 
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r 


34(1 

X7*X1 

44409444 

T4=Y9 

lltldll 

C0T01764 

44449424 

1641 

Tf*T9 ‘ 

miCTii 

X 7* X 1*  (SOOT (N?**2-(N2-Y9)#t2>l-N2 

49449491 

1X9 

_ .. 

▼»tV9 

i mini 

COTOl/o- 

44449494 

1WX4-X9 

tiH«m 

Ti*Y9 

■inui 

114 

C0T0174C 

1644 

TF1 X9. 91. X 1/21  GOTO 1734 

•1199 114 

1994 

14**9 

••••9121 

X7*4 

iimn* 

Y4«Y9 

44449199 

YT5 

umiM 

1730 

*4**9 

44449194 

17  *»r 

•vi«9tn 

V9*V9 

44449 144 

- 

1794 

1 

i 

u 

mm. 

323 

C «£N 

49949244 

c 

*9!ilCTTI77l  TT3  5i2VfT~ 

TI 119 211 

4SSIGN  1774  TO  1SM299 

..  - 

CO  YO  2793 

■ mnm 

1774 

CONTINUE 

49949234 

J7T 

IFCIE.GY.OCOYOIIIO 

iimm 

C 

ASS1C9  1794  TO  SN3299 

4999924 

ISSIW  i/W  to  ii»i2V 

GO  TO  2914 

99999294 

— “■  *'  ' 

1794 

CONTINUE 

•3119771 

333 

1480 

CONTINUE 

44999  241 

C 

ISSICN  1444  TO  SN9144 

•SSIGN  1904  TO  1$M914 

GO  TO  me “ 

11119299 

1404 

CONTINUE 

m — 

C 

TTSlGllTtlf  YO  SWIIl 

nmni 

4SS1GN  1414  TO  ISfctCl 



■ 

CO  TO  9I9C 

•■ram 

1414 

CONTINUE 

XT7V‘ 

uu* i i XUt 

uraiu 

391 

1434 

G0T0219 

44449399 

— - . 

nw-Mnrr  i«9? 

11114391 

1492 

99999379 

6010*1490 

99949399 

1490 

Mi.r  ii62 

99149399 

J<i9 

1462 

7 0NN4T  (IN  ,14X.“**d:77U$I0N~C0U7UNC~ IN-- N0NI20NT4L 

999999|9 

L— sf*ik‘ncK**-) 

99999919 

GOTO!  490 

99949929 

3S  j 


I.T 


vrr 


1992  709N9TUN  ,|&l  •"’DIFFUSION—  COU*.  INC- >IN— VENT  ICAL 

'4— \ 

1090  MINT  272 
MINT  272 

NE40I  ~MTkKtfEU.-,40>t2,N,Cl.T,$.Kl,C7,X5,03,09,09,V9,04,07,09 
eTl^. WfT.tT7W,  69 . oi  ,i9 , 04  767 . oi  • »•;  f9,"i 

163 

AIAO  1C1  •42,09,05.07,01 
MIN?  107  ,*2,09,09,07. 09 
*r&j  ill  ,07,09,09,07,01 
MINT  197  ,07  ,09,09,07,09 
J ,N»  Cl  ,T  ,DW,  fj,  03 
MINT  109  *N,C1,  T,0UN,X5  ,03 

ic*o  i is  .ys.og.s 

MINT  109  ,Y9,09,S 
*f*0  193  ,XL«0,L3,03,L9,04 
MINT  109  .XL, 0,13.03,19,04 

— mjrirs  m , rrn  ,7nw  ,st 

MINT  149  , T9, T9, It, Cl  G2.CJ 
S7-S  ' • 

ICUMC7*  1 

00  2723  UDUN*l,I£U7 
IC«ICUUN-I 

"■TFIIC.W.B71CTT0  T7TI 
X9«X5 
Y1«Y5 

11*13 

L2«L9 

IF (09 «LC .1 £) GOT02990 
GUT020ID ~~ 

2090  Xl*X5M03»tl£-09>> 

17105. LT.lt)C0T022l9 
2040  17(07. IC.1OG0T0219I 
C0T0212C 

2140  Y9«YS»(04VlC«07l> 

"iriOl.CT .irrCQTB77TT 
2120  17  (09.1C.IO  CO  TO  219  9 
G0T02194 

2194  L1*L3M03’U£-09H 

17 (0S.LT.IC) CO TO 2290 
2160  17(07. LC.1CIOOTO2190 
"OTT62299 

2100  L2*L9+ (09* (I C -07  >) 

IFI04.CT.lt) GOTO 2240 
GOTO2270 

2214  X9*XS»(03*05> 

60(0204^ 

27 Jo  Y9«T5*iOl*0fJ 
GuT" 7129 

2298  l l»l  J*  10  J*05  > 

G0T0214Q 


44449934 
•4449994 
•444 49 94 
44499944 
14144974 

*4 144494 

44444949 


44444990 

44444444 

44414414 

44444420 


44494494 
41819494 
44  44  4 4 44 


I HI 49494 
44944744 
14184710 
44494724 
44449734 


44449774 

14489314 

44404421 

14449434 

44449494 
44449444 
fill  W4 
44449444 
44449494 
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ttn  

MW  im.H.WKBTKWI 

I.TMIN 

mi  piimw/ii 

NHW  1 
E*I"T — 


iwitpiwii«"»T"rr 

■•C/M 

— n*r 

wtim ti  t,i 


min  t 

"Tpwin- 


JWWgMtW 

mVMFT 


!WT 


T3*-Y2 

irrn.cr.MunBfui 

60  702930 

tpt 


im 


n»"K 

60T029J1 

irnranfnwfCTTomi 


"TOT 


f2«*t*Tim**2l-ll2-C11**2» 

TJ*-T? r 

triit.tr.shiKiTotMi 


wrr- 


~mr 

2371 

2 Iff 
2990 


MITK21  1 2,Y2 

-prnmrnTrr 

60102370 

rniirr  iztTr 


CONTINUE 
WTTTT7T  rr»6.0 
IFIB9.6T.0160T02660 

nmn 

B7*-IT-U0SIY9/mN<Nllt  11 

THT1  

60T0267I 

i nn — — - • 

17*Cl*MNr  l***2- IY31 **2 » 

rpyi 

CONTtNUE  

1 IS3I61  2611  TO  S*2f|| 

ISSUN  2671  TO  1 9N200 

“¥TTTU  T7W  

i COOT I OOt 

If  HI. 61. 1 110102  fl 


IINMN 
IIM6930 
IIIM4I 
>0006390 
0F5W1TI 
00006900 
• 0116791 
00009000 
•0009910 
00009720 
11119160 
00009090 
00001061 
•00090  70 
00019000 
00009090 

iromn 

00009110 
00019(20 
00009130 
00009  6 

. 00009190 

•nmn 

00009170 
00  009100 
00009190 

oooovno 

0 0 00  9220 
11929271 
00009260 
00009191 
00009260 
00009270 
00009200 


690 


-mr 

2710 


03SI6N  2690  TO  SN  3290 
l»Ilf  2090  TTT  130329 
60  TO  2910 

CONTINUE  ~ 


00009900 
10  009310 
00109320 
90099730 

99999199 

00009360 


01109371 

00009300 

11116398 

00009600 


CONTINUE 

199199  2710  IU  II 
60  TO  3300 

CONTINUE 

mm  271U1A  SN6011. 
0991 WW  2710  TO  I3M60I 


60  19  6191 
2710  CONTINUE 


2710  60T0210 

mi  hiiw  l 


669 


JCHIM0  I ’ 

tININO  6 
DOf«9tlM,7l 
1100111  N,9 


67 r 


nwm  ITT 

«miJI  M.T.S 

02* SONY I 11*171 IT-T01 **2* 

IF  W.70. 0 1 60ro2  Ofl 

60T02000 

>171  B!*!/l 

2000  MNITM61  01,02 


or 


7961  CONTINUE 
C 60  TO  SN2900 

~60T  Iff  Trail, <1770, 1(701 
2910  Of  MINT  3 


NR  61 NO  6 

OININO  1 
97  9160  2 
•( 001JI 
93  2979  I'ON’l.Nl 
0E00IJI 

wmm 


691 


2970  CONTINUE 

C*U  'UTTFUyNITV.Nl.l) 
00  3HII*i, Nl 
TNOTVCI,  (I  ■( 

1011  CONTINUE 

cou  not2iriinoti1ni,i> 
CIU  MT779niN0Ttl,Nl,Nll 
02*6.20310 
tf  *7. 71920 
19*31000 

ini  rwrnirmrrc  ,cn.*wn 

3079  FOONOTII9I1N  ,613.91/11 


1 1 1NI  T0I1  ION  • I COL  I ,t  COL  *1,61  ,tlOV*l,NH 
1 1NOTOdlON,  ICOL  I ,IC0L*1,61 

-TiwnDinnw,  icoc  r ,icot*i  ,4i 


00009610 

00009620 

00009630 

00009660 

10019690 

00119660 

00 119671 

00909600 

00109660 

00009900 

•tWflo 

10009920 
00009930 
00009960 
00019970 
llll 9900 
11119390 
M 019610 
00109610 
II 009621 
11009630 
•I 119660 
Mif969l 
II0I966I 
M|0>670 
1011 9610 

01IV97IO 

101197(0 

10019770 

01009701 

70009791 
01103100 
•0009010 
00009020 
00009030 
00009060 
11097190 
00009060 
00009070 
00009000 
00009060 
00009000 
70 09 fill 
00079920 
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nSUlMfMl  ' n" 

0031SIJ*1,N1 

rrii.EQ.jfeQTOiui 

*l*SQRr<<XNAI0<t,2)-XNATaU,2n**2»(XHATCm  ,3)-XNITOt.J»3n”2) 


6QT031SI 

mi  nuTNfT;ji*AiocTwsn/xiiiTDn#ii) 
ill!  CONTINUE 
mi  CONTINUE 

CALL  NAT iNVt XNATNt XNA  TNt  N1»N1) 


uu 

XN0TIIII *11*1 
“ 3T7I‘CQWTIMF 

00  3121  If *1 «Nl 
00  3&7I  J9*l »N1 

— imiiif  ,rr-TwiTmf.n*xNATN(if,jf> 
3171  CONTINUE 

JI7T  CONTINUE  * 

31 If  F04NAT<613.f/l 


i *«■ 

0032201* IfNl 

iuvwrrrTMmt  ~ 

32?D  CONTINU. 

D032li:*I,Nl 
NEAOUI  01,02 


laimu 

innw 

mhum 

mruii 

IIMMtl 


uiimi 
MMUn 
1 1 nsm 


nrnui 
•i  mi  si 
mitui 


81*UlE-St«X'J)*(M»lf) 

11  *1  lini7-lll  1 

IF  <X.EQ.X9)CQT03ttl 
- " TFTY.T®.TS1C0TOJ72I 

IF<X.lT,X9)G0TO3fll 


M-«TtNt(Y-T9)/fK-XfM 
" T7*-rsnMTJ11 
23*C0S(T3) 
tOTOTTII 

3S20  T 3* AT  AN! IV9-YI /I I-Xf 1 1 


mum 
MIIIIM 
1 1 mm 
•limn 


22«SINIT3) 

— m-Tcosrni  i 

GO  TO  37§G 


IF  fV.IT  • T9) GOT 03 Ml 

” “ ~I7^TSTNTTJ1 1 
23*1 

GIT037IO 

mrz 

23*1 

3 720  IF  IX»GT »I9)GOTQ37t| 


I3UMM 
iiimii 
<iimfi 

• I II  MSI 

iimtti 
iinttri 

inrnn 

Mllttfl 

••••1711 

11011711 

iiiii7*i 

1 1 lit  7 11 

"lilliNI 

il lilt ii 
Mil  iff! 

iiiiiill 

IMIM2I 

iimiii 

iimiti 

iimifi 

urnin 

••HUM 

lllltlll 


23*- 1 

GoTOJ?|0 
3710  22*0 
n*t 

3700  12*11*22 


BS*8S»82 

•6*14*13 

8#  *SQNTI IBS**21* IBt**21 1 
3030  CONTTnuE 

IMOS.EQ.M&OTOJf  70 


IF  IBS.GT  *OlGOTOJf2| 

Tnw.^r.Ticoraifii 

TS*100MTG*S7.2SSI) 

GOIffLIlJ 

3900  TS*U0-<T**f7.2fSI) 

cunyvm 

3920  IF (B6.GT •01GQTO39S0 
ri«36MTIiiS7.2Ull 
GO  TOt  010 

39S0  TS*U«S7.?ff| 

GOTO**  0 10 
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j*rc  xpm.CT.iitftra 


JWM 

omrar 


Mil  T?»M 


ID 


M12  F0RNATI1N  •-NACiniC. 




tn 


ml  -COOR0INA  TC*MU.ll 


Mil  FORNAT  UN  ,-T-COOR0INAT£*-,G13.6l 

pmniTO 


iiiinti 
IINMII 
IIHMM 
I I •••*?• 
Ml  149*1 
iiitMn 

mimi 

miim* 

•0II7III 

IIM7I1I 

11117121 

IMI7I1I 


<i « I37F0RNAT 1 -C00*01  N A TC**,  4 1 3 . 5) 


MM  F0RWATI1N  ,~22-C00R0I MATC-%61 3. 51 

FIWT 


5152  FMNATI1N  , IX,  100*  MCA 

r mom hopt 

MINT  5152 


(DECREES)”) 


IIM7III 

•1112121 

IIM7IM 

' MM2iM 
II l| 2111 
IIM212I 


6012  F0RNATI1N  ,5 I1N  ,613,6)1 

••••7160 

•••17150 

62* 

■HI 
• 1*0 

mma 

••••7171 

■HI 

•3-1 

•5-0 

MM71M 

• IM716I 

• •••72M 

630 

1^ 

C 60  TO  SM61M 

M1I721I 
• ••17221 

40  IQ  I1U6 11 ,11111,27111 
416«  REMIND  2 

•0M723I 

Kk.ni.Ju  - 
REMIND  5 

II VI 7 261 

nr 


f*-i 
om*i>i,«i 
TrTTE.cq.lttOf 06211 
REA0<6)  01,02 
IMBHI  »J 


IIII222I 

IIII22II 

IIM7260 


601052*1 

izii  Wwizr  iitoi iBr.Yj 

5 251 


RT 


F2*S0RTIU**2*02**2) 

T7«icr-Tll 


fl»<l2MU 


F5»iFi-iL?)*mocnri*m/pzm 

F6-IF2-ll*IAl06l  IF2*Ll)/02> ) I 

"F2.isa*T(n»»2ti5f»*2»-rj*iiw6  rrrro«*tir3»»2»02*«2> » /o*»  > » 
FI -t SORT I F6-*2*02**2) -F6* IAt06IIF6*S0RT  4F6**2*02**2) 1/02 II) 

M«iwm?wimw?i 1 — 

FI-SMTtLl**2»0l**2> 


Ql-<F«-Ll*AL06<<Fl*ll)/01)) 

mnwr  in**2*o  -ft*  iaioc  nr  otqrt  if  j**24oi**ti  i /on  1 1 


Q3- (SORT 4F6**2*0l**2)-F6*  tit06HF64SQRT|F6**2*01**2l ) /Oil  t ) 
" A7»rl*M»'rvi  * ~ 


•III231I 

•1117320 

• •••2330 
IIII23M 
•••12350 

nirnii 

IIM232I 

MMFSIO 

MII2MI 
••II2MI 
‘ 11107621 
••••2530 
IIH2550 
••••2554 
0IM254I 

• •••2520 


65*0  CONTINUE 

«mr^* 


••••2550 
• 1112511 
••••2510 
IIII252I 
••••7531 


Ql-0*IXl/R«) 
PRINT  222  ‘ 
PRINT  272 
“FflWT 


•••0255a 

00107561 
•Oil  7570 


5552  FORNAT I IN  *2X  • "TRANSFER^-*-*  -♦-♦-♦-♦-♦-♦-♦-FUNCTION 

MINT  272 
' FIHT  II77  ' 


5572  FORMAT <1N  ,0X,  "TRANSFER  XNOUCTANCC 

r 


TRANSFER  NEST STANCE*) 

MINT  5502 

wrTwiron  jimommiT 


PRINT  5562, XN, 01 
5552  F0RNATI1H 


nwui’r 


12*17620 
•••• 7631 
III!  2461 


PRINT  5602 

6402  FORMAT!  IN  ,"UPIN  CTRCUIT 


PRINT  272 
"FOTOnOIT” 


••••7710 

•0007710 


6616  FORNAT (IN 

— ■ — rrrr 


00  6 721  J0UHNT -1,665 
T7-T7*TI 


nr 


TF  ffTTCT.TO)  W to  6721 

12- I5f  <t-Sl*CXP<-Gl*T7)MG?*CXP<-G2*T7)) 
"T2-t2*t6*lfGl*4li*CXPIt-Gi-G3l*T7)) 

13- 16*  162*631  *CXP<t-G2-G3)*T7) 

~ .5*12-11 ™ - - 

_C7-0(*I6*  !£*PI-G  1*17) -EXP (-62*17  ) )•( 1-CXP 1-63* T7 1 1 

n -wh 


C6-E7-C0 

PRINT  6717  ,7  7 , E 6 


6712  FORNAT I1H  ,613.6, IN  ,613.6) 


672B  WnTTnuc 
5721  CONTINUE 
FRTST  27? 
PRINT  6762 


•0M773I 
MM  7761 
•• 00 7751 
••••7760 
•••07 771 
OOM  7700 
00007761 

r IMITTlO 

00007121 

••••till 

•00I7I6I 

•••07O5O 

rtwmi 

10117000 

•• 01 70*1 
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*4* 


01«J 

”c  w to  miu 

CO  TO  ISMIl,)llllt271l) 


■'mini  Ct€ 

WIT  MJ2 

—\nz  rmm  nr,*nry  rao  srimrur  tots  contain*) 

PNINT  MW 

TWf  nwmintiUB  mm  cmkespond  to  ws-i 

POINT  MW 


r miw  i » in 

POINT  MW 

cm  rnwrnn  v-ro  ie  sow  that  tnct  ik  consistent-! 

POINT  MT2 

wnmTn«  w r»c  cconew  rou  m evaluatinc.-i 


MM79M 

iwwr 


M— If 

mam 


1 

TWC*Jt 

•••••no 

f. 

ill  ntiwue 

IFIMIC.E2.J1I  CO  TO  211 

•■tiff  21 
••Horn 

\ 

00  741  JZ«1,NC 
TENP«A  IJ1.J2I 

ffifirci 

•III02M 

CO  TO  MB 

Tirwi 

NC*K 

oo  n jui.wo  — • 

21 


TNP1M. 

* -00  in  J7«J1.N0  " ■ ■ 

TNP2«CAIS)A<J2,J1)) 

— “TWPzoisurJT.jni  ' * 

IflTNPI-TNPlI  121,121,1211 


IIIIH2I 

nirnn 


•••Mill 
IlflUEl 
• III  MM 


lAOEMJlI'lAOEKMIC) 


211  TENP-AUltJH 


00  221  J2*l, NC 
til  MJt,J2TMTJl,-J2T/TtOP- 
.....  00  211  J2«I,N0 

1FIJ2.E0.J1I  CO  TO  111 
TfNP.lfJi.JU 

Al J2t Jll •• • 

W f¥i  JJ»i,NC 

241  A( J2,J3I  «A<J2, J3l*TEN^*AIJ|» J3I 


COflTINUC 

211  CONTINUE 

itnii'ini 

00  301  J1>1,N1 

“ "win  J2*ji,«r 

IT (LA0EL I J2)  *NE. Jl)  CO  TO  J2l 


IF 

00  TO  |»| 

Tmwnm* 

311  00  301  J3«l,Nt 
IU),J1I«I|J],JZI 


MIIMM 
ill I lift 
ffffiO?l 

•Ill  MM 

liofiiM 

•I  MUM 


IINH1I 

um 

iiiiiMi 

iiHiili 


::ui?u 

•mini 

iiiitiii 


Figure  58,  DIFFUSION  Program  Listing  (Sheet  8 of  9) 
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t% 

•• 

Ml  *tJ),J{C«TCO 

LIBEL  IJZMMBCU  Jit 
Tfl  C BIT  I HUE 
M 10  000 

* WIW  «flTTI«M 

00  Bll  !2*1»L 

ninm 

IINNM 

0000*091 

•••BOOM 

HiMin 

I0009BII 

■ - 

VII  EIII.Ill.ini, !?l 
Ml  RCftMM 
“ ' CTO  

mm 

••MUM 

ihhiii 

s 

SUBROUTINE  NRTINVM,B*IROM,ICOLt 
■*  OTHFHSTOK  HI  RON.  I COO ,Bf  IRON,ICOL» .CU.M 
00  HM  t«>l  •IRON 

~ 00  T5T»  , I COL 

•*(!*•.*» 

00000200 

BBHOZtB 

BBBBBZZB 

0000I2M 

BBBBBZBB 

11 

J!JI  C01I IMA 
CONTINUE 

CHI  UTmtl.B.A.C.  IRON, ICOL. IRON, ICOL.ICOLI 

RETURN 

CKO 

mum 

0000IZM 

•IIIIZTl 

BBBBBZBB 

BIMBZBB 

s 

SUBROUTINE  NOT 2ERI XM, 1R0N,IC0L  1 
OIRFWSTON  IN  t IRON, ICOLl 
00  MM  » IRON 

“ 00  Jiff  JB-l.ICOL 

XNtM,J9>‘0 

•BBBB1IB 
BBIBBIZI 
IIIIIIJO 
MM  a Ml 
IIMI19I 

ini)  LUli  mm 
1211  CONTINUE 
RETURN 
Eno 

mi  mo 

00000100 

BIIBB1M 

Figure  58.  DIFFUSION  Program  Listing  (Sheet  9 of  9) 


150 


pkocram  aptun!  Input,  output) 

c APERTURE*— -A  PROGRAM  THAT  CALCULATES  THE  MAGNETIC  fnU  THAT  ~~ 
C PASSEi  THROUGH  AN  APERTURE*  PA  FISHER  ILK  f-tll 
C itNtML  ELECTRIC  COWPANT  tit  KOOKMN  Iff  PITTSPIElt  »•*«*  tlttt 

c phone  mu-m-nai 

C QtrtJLOffO  UNDER  CONTRACT  PIHttWIHI  -NSRP^LttWT  OTNRHte*  t 
C THE  PROGRAM  NLA  OS  DATA  FA  ON  AN  EXTERNAL  FILE*  THE  MANE  IT  WHICH 
C WILL  AE  KE  QUEST  EO  WRING  UftCtfttON,  ~THE  INPUT  NAT  A Pttf  SHOULO 
C WE  CONSTNUCUO  AS  FOLLOWS* 


LlMt  MiNlEA  ll  KAfTA 1 2A 

tt  LIlLt't ANtW  -- 
II  NEXTtAN6N 

At  oxttfr 

«•  oi 

•i  gpAitPitfp; 

71  TPA,YPI,TPC 


C 

c 

~c  • 

c 
-c 
c 
e 
c 

r 

O' 

c 
c 
c 

• *t*  we  nunqcas  i 
c 

o-  «*rr*,t*  tit  **t - coon ot nates  m-Netfuruf  tnc  center- cp  twe 

C APtRTuAE.  IT  IS  L0CATE0  IN  A PL  AIR  PAftALLAL  TO  TNf  XT  PLANE 

i u ANJ  L2  an-  THE  LENGTHS  in  NETENS  op  INE  AXES  of  THE  ELLIPTICAL 

• ♦ it  AT  L*  tl-NAJOR  ANtS  ANfr  ■ tfWWHt  Atttl 

C Amah  is  THE  ANGLE  that  THE  MAJOR  AXIS  F THE  APE A TUNE  HAKES  NITN 

- -e  -INC  A AXIS.  • OtONttS  IT  P ANNUEL  -tWE- 


II  OA 

iM  Tt  - - 

Ul  PXl,PTI,P2l,PX2,PT2,P22 
in  FT)iff»iF»<mtmiFh 


»Nc 

C 

NITN  NCSPfOT  TO  M-t  AmTi  I HtNEtf  PANALLCt  f 1 


ES-THE A£  IS  A REFLECTING  SURFACE  PARALLEL  TO  THE  APCNTUNE. 


-fr-  Jf  2 eOONOlHAie  tf  THE  REFLECTING  HHFWtf — ENTEA  iUNNT  tALUE  IF 

C 01-0. 


IIMIIII 

tint  tot 

mum 


C 03»1»T lS-CAlCULATc  THE  FIELOS  OVEN  A PNESCAIIEO  VOLUME  I NS IOC. 

-c- 

c 


•1111*11 


C 2PI-2  COORDINATE  AT  WHICH  CALCULATION  SHOULO  ENO 

I DW  iHONLHcNT  UK 

C TPA,TPIfTPCfXPAfXP|,XPC  AAE  SlNlLAN  F ON  X AND  T COORDINATES 
-O-HtfW-  OWHHT  WAVES  IP  Oil 

C 

-O-PhM-T  AOUi-ATc-  flfil1  IN  OPWCNIOAL  OAAAWWATHi 

C 0*»-C*T  A BULAT  £ IN  RECTANGULAR  COORDINATES. 


05-1-Yc.S-CAlCULATl  THE  FLUX  LINKING  A LOOP 


inutr# 

ilium 

•••Hill 
tttllltt 
••••1111 
iHtiiet 
mill  si 


-i 


fXl.Ptl,—  -PYA,P2A  ARE  THE  COORDINATES  OP  FOUN  JOINTS  THAT 
0CF1NT  THE  cOOP.  THtTTIUST  CO  RROUNP  TWT  LOOP  IN  CO  NS  0 CUT  I /E 
OhOlH.  ADDITIONAL  LOOPS  HAY  IE  OEFIHEO  IT  AOOITIONAl  OATA  IN 
THE  SANt  FORTAT.  OUHHT  VALUES  ARE  NO?  REQ9IRC0  IT  Olw| 

U1HEHII0H  HH(tr,  tST  “ ' ' - - 

01 HlHSI OH  T I6A (IS) 

DIMENSION  PATH  At!  JT*- 

<*  AL  Lld2.HUl.NU2f  NU3,MU,NUX,NUV,NJ2 

PRINT  111 

C«-tkI*bE  CONTROL  FORNAT  STATEMENTS 

no  fornatuh*)  “ * 

Ill  FORMAT  I INI) 

120  FORNAT (IN  ) * * 

122  FORMAT <1H4> 

T2S  FORMAT (IHft- 
OUirUI  OATA  FORMATS 

1SI  FORMAT  l*£l*.  Sr 

DATA  HE AO I NO  FORMATS 

1M  FORMAT (-  APERTURE  CO OROtNATES—lx*, ItlT.Ii"  METERS*) 

1<*I  FORMAT !*  T--,1E12.J,*  *1£TERS*I 

1IJ  FORMAT  ( " fx*,ltl2.3»*  CTERS1 

111  FORMAT  ( " APERTURE  01  HEMS  IONS— MAJOR  AXIS*-,  tl  12#  I, 

(*  METERS*!  — 

Ill  FORMAT  I-  MINOR  AXIS--,  II 12.  I* 

r METERS"! 

Ill  FORMAT ("  APIRTURE  INCL IMEO*.  1E12 . S. - DECREES  FROM  X AXIS*! 
ire  FORMATi-  tXTLRMAL  MAGNET  SC  PtttO**,  1C12.I, 

A-  AMPERES  PER  METER"! 

171  -FORMAT  1-  AMO  WOtUNTintW ttlNggS  FUR  fXf-t  )»Ht 

u;  furmat (-  there  is  ho  reflecting  smface-i 

111  FORMAT  (-  TPtRt  I**  REFLECTING  SURFACE  LOCATEO  AT  2«“, 
IU12.3,  " HcTtRS"! 

Ill  FORMAT!-  LOOP  MJMttR  *.II) 

1)0  FORMAT!-  LOOP  ARIA«-,1E12.S.-  SQUIRE  METERS") 

tl2-FORHAi  <—  TOTAtrPLVir**rltTr.rf'“  - WCURI-T  - - 

1)1  FORMAT!-  OUT  OF  OATA-) 

211  FORMAT!"  POINT  X V 2") 

>24  FORMAT! I, ,JE12.J! 

Rt  AO! INF IlE • 23t, tMO- tl*l)Ll ME ,XA  f VA,  2 A 
RtAU  <lMflLE»2tw,EH0sltfetlL  INC, L 1 , l2  , ANGA 
Rt  AO  lit  #L,-A#YAyr«~  


imisff 


mmtt 


I lift  III 
I JII1I2I 


IIII17II 

mum 

IIII172I 

mum 


nwtm 


•iitirtr 

IIII1I2I 


mtiiu 


•••um 

muni 

mum 

ilium 

ilium 

•limit 

11112111 
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w; 

1*3  FORMAT (fctl*. II 

RIAO  113  ,Lt,tT,Rl«r  - * 

PRINT  10I,)IA 

'ARINA  199, TA - 

PRINT  151 tZA 

>31021 21 
■Will  3 > 
MOIIIAI 

PRINT  100  ,L  2 

PRINT  149, AN* A 

IIIttOTt 

C RE  A0(INFRE,  200, ENO«lf*d>UNt  ,*€*?, AN«M 

RtAD  113  ,*€*t,AN*N 

MH2>» 

IK 

PRINT  *115 

"PRIMI  irffHtFT 
PRINT  179.ANGN 

MM211I 

PRINT  Ilf  - 

timm 

C A A(K!Nf  IL£,2M,ENO«1900>lINE,01,Q2 

>•992131 

115 

Ac  AO  U3  ,01,02  

- 

ft  9 IF(Dl)  2*0,?M,m 

•VI 1219V 

2Av  AN1MT  11) 

>>»2tf> 

ARInT  115 

OHIZtlV 

215  G0T0255 

it: 

f*R  Allot  l#f,02 — * 

tltttlW 

AnINT  115 

>>»215> 

295  CONTINUE 

*>nm» 

C R.AOIINF  iLt  ,2M,£NO*19M)LINE,03 

>>»221> 

C RtA0UNfKl,2ll,€N0»l99IILlNt,rAA,2AI,2AC 

>>>>2226 

iZi 

C K- AO  ( iNf  KE , 201  ,EN0a19ll)  LlNi  *YAR*TR>*YRC 

• •>>22  3> 

t RtADtlNAitf  ,2t*it*R*t9RttttNt  ♦lAAfRAf,  MAC 

DfMZtNt 

C Rc  AlMlNFRE  , 20l,£MO»19fc0)LlME,OO 

Ri AO  lv3  ,03  — 

R-AU  133  ,ZAA,ZAt,ZAC 

DM2230 

IK 

NiAO  Id  3 ,TAA,TAtiTAC 
R AO  103*,AAA , A At, ■ AC 
RiAO  11 3 ,0%  

PI«J.l*i5l2o5 

• 3M226I 

C»U  SNAAt(ll,L2,An,tn) 

0 1902271 

1 J5 

IF  (Oil  1952,  1952,2211 

DM2200 

22)  J If  (0*.l223l,229l»Wt9 

> t>V229> 

22.1  AKINT  2232 

(Zti  FORMA H**  * — -Y  I 2 N-* 

wttntt 

; a n-y  N-z-i 

>M>232> 

11. 

GOTO  22(c» 

••>>2331 

2..t>  P*INT  22)7 

> 310234  0 

ttjf  FORMAT t-  A - Y 2 N TOT 

»M235> 

4 LAT  LONG**) 

IIM23M 

. . 

22 /A  PM  NT  1 2d 

1 MM3  At 

1-v 

Goto  2453 

> QM23I0 

22.4  P-JNT  123 

D>>7  3 5> 

2450  CONTI NU£ 

>•>>24  >0 

J*«lf  1«  ( (ZA8-2AAt/ZACIil 

MM2%1> 

J2*IFIA(  (YP3-YAAI/YPCH1 

MM2420 

if* 

- IHH*A8‘MRi»A*RD»*l 

tttttiJt 

□01953  IIU,  J3,l 

DM2440 

001990  I2«l**ftl 

>d>>269> 

30 1953  Il«l,Jl,l 

1 00>24>> 

«Ai»«#AitIf-l)**AC 

>>»»%P> 

Ijj 

YM«TPA*  U2-1IMPC 

• >M24>0 

— — 

-2Al*ZAAi(U-i)ifAe- 

MM  AMO 

1332  SuNTINUb 

DM2510 

*1*9  C6NTIMU* 

• >M2f  1> 

1 753  CALL  MA^fLOIAlGA.ANGN, API, YA1,ZS1,AA,VA,ZA, NEAT, A11,A22, 

• M>292> 

10- 

i HlH(l,MAfl,HA2l,tl»C2i01»02l 

>5MtfM 

3173  IF<Oi)l21J,121d,A090 
I2t1  PRIMf  l2MtrAAtrtAi,2At^MAAlrHAYl,MA*l 

• 1M2541 
•9M2999 

GOTO  1353 

• M9256G 

l^d  -FvP*AT«4M2. 3) 

tttttfPO 

145 

4 1,  . J j:  SORT  (HP«l*MPAliHAZl*NPZl) 

t JM2SM 

IK»9S<*»Tll-*»3<*t) 

>3>tt99> 

*:in  ANGl*9rf-57.2957795MATAN(NPYl/0l  1 

MM2M0 

m-tvk  w»T0  4012 

»>mi> 

4 J 1 (t  A««C  1 • 57  »295  77  95*(AT  AN( 0/MAY  1 1 1 

>j>>2>2> 

If, 

If  t A85(9)**i  1 -ARSTMA21I  1 lull  ,il  1 4 ,40  2t 

MM2439 

i ili  AsG2»9J-i7.  2957795*(ATAN(NPAl/HPZl)  1 

DM2441 

4,16  GOTO  4*  3 3 

»»2»5> 

4 J2  S ANG2>57,  ? 95  7 7 9 5’ ( A V AN<  NPZ 1/ HPAD  | 

03002M0 

4 > 33  I f | hPY  i ) i;5<o4d50  ,iiiO 

»>02(7( 

If'. 

4.*:  SOTO  4113 

01002000 

mr,i  A»JG1«1D-ANG1 

•3002090 

*i;d  scntinu; 

DI027I0 

>.123  IMNPAllildw, H30, 413. 

039327 10 

413)  IF  MPZ  1 1 41o0  , 4 K J , 4 14u 

03002720 

18 

•1«,  ANG2>ANG2 

•3002730 

415  3 GOTO  <>215 

03132740 

4 I-,*.’  ASG2B~AN*2 

>3912750 

4 1 7C  SUTO  4213 

• IM2700 

418a  If (NPZll421w,4l9a,4l9; 

•9002774 

155 

4193  AMG2«18:-ANG2 

DM27M 

i2, J jOTO  4215 

DM2790 

421J  A'<G2"*(  1A9-ANGZ1 

•IQI2V00 

4217  SC  NT  I NOE 

DM2D0 

H«T*S0RY(NPI1»NPA1*NAY1*)*Y1*I#21*MPZII 

D0>2>2> 

IK 

MR  PM  NT  1223,  API, YP1, ZPl , HPT,  ANG1,ANG2 

•01020 30 

1)53  CONTINUE 

0>0>2>4> 

1952  •»-  INT  in 

MM2990 

S R- AO (INF ICt ,2ud,-NU«l953llIHE,05 

R-  Ad  ns  ,05 

• QII2I60 

1 95 

2KJ  If  (05)  14)0,  1400,1953 

>mnr> 

1)77  CONTINUE 

>>> )2>>0 

J*»L 

>•••2190 

1)^7  CONTINUE 

M»29>> 

3 N-A0(lNFn£,23d,EN0*19ll)LlNL,PAl,PYl,PZl,PA2,PY2,AZ2 

DM2910 

Figure  59.  Computer  Program  APERTURE  (Sheet  2 of  7) 
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MM 


1 

\ 

} 

V 


l 


r 

E 


m 


^lOimriLEtZiiteNOMHiiLiNCtnstHSi^iimiPVMn 


IF  INlMMIIMIlMMMlilIN 

C THESE  ARE  INC  SIOES  Of  THE  QUADRILATERAL 
lilt  CONTI HUE 

«>f  «W«t 

iniMin-mi 

Xl4aPXl-PR4 

IJWTI-W 

T«|MV«^V9 

T14"PT1-PT%  

221*P22-P2  1 
242»P2S-P22 
24I>R2«-R2I 

C THIS  IS  A 01  AM  HAL  Of  T Hi  QUADRILATERAL 


231*Pfl-P2t 

Tn*soRT(Ui*»i*rn*m  *m*mi 


TO*$M(  IMI«US«T<»I*TkMAI*lAI) 
1 1 A*ttW  Wtmtt  4*  214» 


ESS 


>ilaSQRT(XJl*X31*V31*V3l*2Jl«231> 

• 

Ai«SQftT<SlMSl-mi*CSI-T32>MSl-T3in 

»MT«3«UMT3ll/2 

Amour  or  1 3i>mi  » tst*  mmsE-isi  > 1 

AREA»A1*A2 

C 1«m  ARC  INC  N10P0IHT3  Of'TNC  CMOS  Of  THE  QUAORI LAT ERAL 
KfHl«fU»Ml/2 

• • 1 • VfRt#Tt*Tf  I ft  ' — . — — 

2PNlaP2 1*22 1/2 

KWI*m-WTTl 

V2N2<fYA-ni/2 

z*m*n%-i%srt 

XPN2l*XPH2-XPNl 


29# 


in 


2PN2l«2M2-2fNl 

T PA**  SQR  T fxfw  rmPHTl  ♦ YPH2t»  T PH2 1 ♦ 2PH  2 1 • 2PH2 1 1 
C TRCSt  ARC  TH£  :ONP0HEHTS  OP  THE  HORN A L /ECTOR 
nui»m*»2*Tjmn 
HU2**I2l*£32tf32*22l 
“NU3al2l*Y32*XJf*Tl  t 
T.NPaNUl*NUl*NU2*NU2«NUSsNU3 
tWhSQRT  ITCNP1 - 

C THs.Sc  ARE  THE  00  HP  OHO*  . Of  THE  UNIT  NORNAL  /ESTOR 

— • ~NUX*NU1/NU — 

NUY-HU2/NU 

HefHHt/Htf 

3 SSI  CONTINUE 
- -JWt-W  lltl  N2-i,t*,t 
3S7I  003110  Hl*l » 13*1 

W x 1 * X 21  ♦ tit 

Y°3*PT 1*121  * CN1*1I  / 12 
— 

XPb*PX2«x32MN2*l) /12 
TM*PYf«fH*WMWtf 
2PA*P22*232«IN2-1»/12 
- - -»Pf*PX4-XA3«t«l'-ll/t2 

YPf«PT*»-f  A4#  IH1-1I/12 

2PMP2*-?*l*fHfi»/t> 

XP|aPxl-X14*(N2-ll/12 

— • - YPf  PTl-fl4MN2-l>/l*  - 

2PA«R21-21A>IN2-1»/12 

— - - — - - 

Y79»YP7-YP9 


I«MXP|>XP| 
-f«4*YP4-YP4  — - 

294-2P9-2P4 


XPL-XPI- 
YPi*YP| ■ 


bMJHI /12 


2Pl*2PA-294* (Hl-ll  /12 

4 tNPYl»HP21fLl|L2f01|92l  * ' 

- MHPNI»Wl*HUW«HPY|*HHYMtP2fHHl 
TNESc  ARc  TH-1  NH*S  AT  THE  VARIOUS  POINTS  OF  THE  3U AORlLATERAt 


C lHc$c  ARC  Thc  OISTAHCcS  TOP  10  IOTTOH  AlONC  THE  QUAORILATERAL 

- T94*«N«F*tt4 

4IIQ  CONTINUE 

Slid  0039SJ  N2*l|  1 5» 1 ' 

- orLTAi«Tvu<*2>rir 

C THIS  ‘.YALUATiS  FLUX  ALONC  THE  LlNcS  1H  THE  OIRECTION  Pi— «P2 
« ANO  PW-fPJ  — • 

PA  THaNHI ltH2) »HN(2,H2I •S*HN < 3»H2J *HN ( A , H2I •»♦ NN( 9.N2I ♦NN(4, 

CN2) M*  XXr7gNmT*in(tpmM«lll(9lIIII*HN(ll,<l2)*l» 

& ilH( Ilf H2)aMH(  12(N2l *9*HN(  1 3(N2) 

PITH  A (N2T*0.  3*  BELT  A1*PATH 
3«0  CONTINUE 

0:iTA2«TPH/l2 

C THIS  .VALUATES  TNc  RESULTANT  FLUX  IN'  THE  0IRECII0H  Pi— spa 

• x-nwr  ff-^PT  — - ’ - 

HT  OT » PA  THAI  1) ♦PAT HA ( 2) *5*PA  THA< SI  *PAT  NA  (A)*  l*P  AT  HA  I SI ♦ 

^ PATHA  HI  *S  ♦PITHITTI  *2YPATHA  fl)*9*PAT  HA(S)*PATHA  ( 10 1*4* 


•MI292Q 


1 1112931 
IHttlAt 
49112999 
HHffli 
•NI2ITI 
•iimii 
1)112991 
•totstto 
tlllllll 
ItHSt  (I 
49493434 
•I ««3« A I 
9)M3|S9 
•9999949 
99993979 
99999999 
999999  99 
99993199 
99993119 
99993129 
99993139 
99993149 
99993199 
99999199 
9)993179 
99993199 
9)993199 
99993299 
93993219 
99999229 
99993239 
99993249 
919932S9 
9)993249 
93993279 
99993299 
999992S9 
90993399 
99993319 
99993329 
99993339 
99993949 
'999939  S9 
90993341 

99993379 
99993399 
99993399 
99993499 
93993419 
99993429 
*9993439 
9)99  3449 
9)993499 
99993499 
93993479 
9)999499 
99993499 
99999999 
93993919 
•9999999 
10913539 
99993949 
99993999 
9999)949 
99993979 
9)999999 
99993999 
99999999 
99993419 
9 99934  29 
99913439 
99991449 
99993499 
99999449 
99993479 
99993499 
9)993490 
99999999 
99993710 
9999)729 


99993799 

9)933740 

9)993779 

99993799 

1)143799 

93993919 

99993929 

9)999939 

9)993199 

99193949 

99993979 


i 

i 

i 

5 

t 


t 

I 


j 

i 

i 
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sac 


3j5 


Jl. 


SIS 


42 . 


1 PAfH#(Ui*P«lH«tl2)»5*PAfMA<L3> 
NT  Ot  • 1 . 1*31  LT6  **l»T  OT 
*1  Ql*Nl  0T#«*PI*1E»7 

• * 

■•415  PAINT  144, JI 
PRINT  1*3 
P*<INT  211 
P<1NT  1*3 
J*1 

pni**  *M,j,**n*Tt,*rr- 

J‘2 

p<inr  Hi,  •,rx2t+rit+2t 

J-3 

PRINT  ZZ3»J,P«S,PTS  ,PZ3 
J*N 

P*lNf  Z?«,J,*«vv*T*T*f* 

P*  1 mI  1*4 

N 3*0  PkINI  Hi, ARE# 
faiiQ  PilNT  152,9101 
*3*3  PRINT  lia 
»UT0  1 55  7 
■ 3o9  PRINT  113 
15*0  P-1NT  195 

i*;c  stop 

LNU 


iUflRJUllMc  Sh#  Pt<Ll  ,L2,«U,  A**) 

R;#L  tl.W 
P1»1.1»1>9Z** 

£ l* 1* (L2/L1 I ••* 

£<r>SQRT(U> 

t:np*lz/li  ■ ■ — ■ — • 

If  ( TEHP.cT  . C . 0 1 GO  TO  JlJO 
IF  (TENP.CT, 1 . IT  50  TO  1161 
: fMMMflMfffAAIAAffAAMAflAff 

c cr«.m,EZ»  a no  cELir?,c2i  arc  itath  lhr#«t  routines  twt 

C . V AlOA  1 £ THc  ELLIPTIC  1 NTfcGRfiS  Of  INE  FIRST  9 NO  SECONO  KlNUS 

— cali  ccuiTi ,tmtio  — - ---  

tJ»l.«J-E2*£2 

GALL  CcLZCTZ « E 7, 1 • | « 1 3 , 1 EWT 

C IMfMIAffMffMAAffffAAfMffM 

com«**Pi*ui/n**3/j 
All«CONl*cl/ ITl-TZ) 

Az?*cdNi»ti  •ii-nTT<T7M  i-tnnn 

A J 3 *C ONI*  (1  fccl  l/T? 

RT  TURN 

3 A 1C  PtlNl  Jlta 

31*0  FORNATf  LZ/L1  TS  NEGATIVE.  TM3  TS  in  £*R0R“> 

Hi  Turn 

3 16 j- print  3iM  — - - 

3173  FORH9TC-  LZ  IS  LAK&iR  TH#N  LI.  THIS  IS  #N  ERROR") 

STOP 

t*0 


a m mi  a 
a mm* 
mmii 

■1113*11 


iiiimi 

tut  sift 


tlffJlW 

1111*111 

timtu 

•111*121 


ttmift 

MNMM 

tmtart 

1111*111 

1111*191 

1111*111 

vnmii 


nutisi 

iimni 

mmii 

11119171 


nmirt 

1111*211 

9999*219 

1111*221 


•m*rct 
• Ml  *2  7| 
1111**11 
1111*251 
ivmivt 
•311*311 

9999*979* 

1111*3)1 

1111*1*1 


•199*3*1 

•111*371 


15 


t. 


2-> 


37 


4. 


SuBkOUTIRe.  NA&FLUlANGf  ,#N5H,*Pl,  TPi,*P  ZA,NE*T,»  U, 

1 L?Z,HPA1,HPT1 ,H*7 1,L1,L2,0I ,02) 

R.fL  *1,<2.*5,<*,lI,lZ 
N*0 

R4U»37.Z95779? 

* Pi  • 3*  1 n159?*3 
A1»o0SJA9UA/R90) 

N7*SIN(#5GA/R90| 

C CAcLUt  • T ION  OF  SHlFlcO  COCRulHAUS  OF  POINT  UNDER  INVESTIGATION 
AP2»«Pl*Rl*f Pl**2 
fPZ*-APl*R^*TPl*m 

tp*«2Pr 

HP*Z»L 

HPT*»8 

HP22«m 

/AA*2# 

C ^AlCJ*.#TI0N  Of  01  stances  *f ROH  aperture  TO  POINT  UNDER  SlUOT 
5133  Xl*ITP*-IA 
51N4  ft*TP2-Tl 
9150  Z-»ZP?-!AA 

51bw  Cl«*AC#AG*TC#TC*2C*2C*Ll#Ll/* 

CZ«AC*AC*TC*TC*tC*rtn2*L2/* 

; Gf  LL  Jl  A T 109  OF  H F ItLO  P#R#LLcl  TO  #M£S  Of  #PtRTjRE 
CUN1>1,*51J3-  - 
CUN*>  1*  A J*5  96 
CON|*2,  b 3*727 

C SAlCJw #T I ON  OF  field  PAR# Lit  u 10  AMES  OF  6PER1URE 
#N#H* (6NGH-AN06I/RA0 
H4AU>Ht  XT'COSl AHAH  > 

H«IN«HHT*3INUHAH» 

v GAlvJl#!  IUN  OF  R0T6TED  C0HP0N£NTS  OF  NA5NEUC  7IEL0 
Nx  «•#! l*HHf 4/ l**PI*Llf 
No* ■ AZZ*HHlN/  («**PI*LZT 
Li  * *C*c 1 

9 1*5  If  <6di(C3) -l£-5»  5A6J, 5174,5171 
91 7y  F i»C  J#C  3 

F 7 « 1 *Z*F  1*5*  F 3*F3*7*FJ#F  4*F  4 
*«T0  5191 
9161  f/»l 


1111*311 

•199*359 

1911**11 
1111***1 
•m**  si 


•919**51 

1)11**M 

1399**71 


•119*911 
• 111*921 
1999*991 
OM*5*l 
1191*991 


1199*971 
1 Jll*91 1 
• 199*999 


1111*1 11 
1111***1 
1999**91 


9999**91 
1)11*711 
1119*7 If 
•111*721 
1191*719 

• 191*791 
1111*7*1 
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•I  41*1  M«C1/CI 

9199  IFi»*S<CS»-lE-9»  9211,9211,92*1 

■ -+f*r  Off>*C9  ‘ 

Fl«C9*ll*C0Nl*C7*C0N2#07*C7*e0N3*C7*C7*C7l» 


• -••  GOTO  un 

99  9211  Fl*C9 

9tH  -t%*TC*Lt  • ~ 

9^31  IF<A9S<C*>-lE-9>  929C,92*0,92*i 

9tM  >9»e9*t9 ■“• 

F»«I*2»F9«9*F9*F9»2*F9*F9*F9 

91  -1W9  9299  " 

9299  F|«t 

9m  C»«C9FC2  - - 

9 If  I IF(A0$(C9)-lt-9)  931W92II.92N 

9 mi  eM»»e* — 

99  F2«C**  I l*CONl*CI*eON2*CI*CO  ♦CON3*CO»CO*C«1 

— - -49M9J2I  

9309  ^*C9 
9121  CONTINUE 


ti*3*K9*«CACr~tT9 

§t  62fJ*Kfc*lC/C2»M.9 

— - - G9*«9»ll/Cr*M.9- 

G**3*«9*VC/C1”1.9 

-tOifW^FKMi.f 

6MM*L2/C2”1.9 

-it «f*WK^l*nrt 

il>]*m2C/C2*«l.9 

• - - mht»Ci*Fi*fc2»rf-oj*fT 

HFTit»*Fl»69*F2-69*F| 

tt - 

fC  9911  IFIOU  97*0,9791,9912 

....  W||  h rtf**tfn+x  • 

NPY2»NPY2««Pf 

“ NHf«H22*HFf 

9991  IF  (MS(NFK)  • • A09*A  OS<NPX2) ) 97  J9  , 9711  »9f  21 


9711  IF  IAIS<NP2>-.0I9«AIS<HP2?M  9771  *9770, 9721 

im  t*  twin  9739  .iritttm 

9731  N*N»l 



•0  GOT 09131 

9 791  HPNf— PWfNPN 

NPY2»NPY2»NPT 



9771  CONTINUE 

*—•9 - fr-C*t€*t-*?ION  OF  C9NP9NCNT9  Of  MAGNETIC  FIEtO  OOTITEO  §AC< 

C TO  TNc  REFERENCE  AXIS 


HPT1»HPX2*<2»NPV?*<1 


- *Pfl«NP22 

90  RcTURN 

• ENO 

0 

c 


191 


110 


119 


• TTtr 


12& 


ist 


139 


C SUBROUTINE  CU1 

9~  - 

C PUNPOSt 

■O  HtOMH  HIFitti  CttlfTIO  HHWH  IF  HW  <INO 

C 

e - USAGE 

C Call  CEL1<RES,a<, IER) 


0ESCRIPT I ON  OF  PARANiTERS 

A < - N00UIUS  I INPUT  I 

ERROR  COOt  khC RC 
IER»0  NO  ERROR 

- K*»l — A*  Nftf  IN  RANGE  -1  TO  ♦! 


FOR  •<•♦!, »t  TNl  RESULT  IS  SET  T9  1.E79. 

FOR  NOOUiOO-A*  ONO  MMK€*€*TARV  MOOULJS  0<, 
EQUATION  A<*AK»CK*C<»1.0  IS  USlO. 

A<  NUST  0&  IN  TNE  RANGE  -l  TO  *1 


SUQROJT I NiS  AN0  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 


METHOD  ~ 

DEFINITION 

mirmi.TiiTTpi«Lri/SQjrnu*T*T»*UACC<*T»**2n , sjnhep 

OVER  T f RON  0 TO  INFINITY). 

EQUIVALENT  ARC  THE  DEFINITIONS 

Cti.  1<A<MINTEGRAL<  1/ ( COSI T I S3RT < 1* <C< *T ANI T 1 1 ••?)), SUHi£0 
tTTTRT  rROW  9 TTT ri/21, 

CClKARMINT  cGRAl  ll/SQPT(l»IA<*SIN<Tll**2l  ,SUNNEJ  OVER  T 

FROM  3 T0PI77I,  MRE1E  <*S0RT  1 1 • •C<*C<),  • 

evaluation 

LANOEWS-  TRANSFORMATION  IS  0SE3  FON  CALCULATION, 

RcFtRtNCE  

~ ‘N.-murNSCNT'-NwriRiciL  cn.  cm  it  iow  of  elliptic  integral! 

ANJ  ELLIPTIC  FUNCTIONS",  NANJOOOR  SERIES  SPECIAL  FUNCTIONS, 
NUNENISCNE  NATHENRTIK  f X.  F,  1919,  PP.  7|-9I. 


MM977I 

0) 009711 

1199*791 

1) 11*121 
• HIM  31 

1999*9  99 

1)11*171 

1)01*111 

919f«ftt 

1)11*921 

1)11*931 


1111*971 

•911*911 

•9999991 
101 191 II 
1)119121 
1)1191*1 
•11951 *• 

99909991 
0011*171 
III09II 0 
MII9I9Q 
19119  111 
I JII91 II 
11999921 
111)9131 
111191*1 
IIII9190 
1099*191 
II0I917I 
•9999111 

0) 119191 
I9II92II 
09109211 
14019221 
09009230 
1*9992*1 
09109290 
11999290 
IIII927I 
•9999219 

1) 119291 
OftltttO 
.CEL  10021 
CE11993I 
CEL1II*0 
CEL1II99 
CEL10 190 

CEL  II 0 00 
CEL10090 

CElIO 1 0 0 
CEL10110 
CEL101 20 
“Of  Li  91 30 
CELIO 1*0 
GEi.10190 

CELIO 1 60 
CH10170 
CELIO 1 00 
Cftift  99 
CEL10200 
OEL1I21 0 
CEL  10220 
CEL1I230 
CEL102L0 
CEL  11291 
CEL1I290 
CELt 1270 

CCUfZM 
CEL10290 
CEL  1030  0 
CEL  11310 
CELIO  320 
CELIO  330 

cnm*o 

CEL10390 
CEL1I360 
CEL10370 
CEL1030  0 
CELIO  390 

cctu*oo 

CELIO* 10 
CCL10420 
CCL1I*30 
>CEL10**0 
CEL10450 
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SUBROUTINE  CUjUUSJU'XOU 
X«.MI  ^ 

TEST  no out us  

6‘i.on. -**•*« 

tFfftC  011,2,3 • - 

. * 

RETURN 

nt  wmt  nw  ■ortw 


t res«i.  crs  - 

RlTURM 

t-fRO-HRTttCUf 

MX*li 


USTiMIXM.t^ 
RRi»ftio»ARr 


- ' ffSt'  'If  ICGURHf 


S «0a  SORT  I AAR  I*  If  01 



60  TO  I 

..  _____  J . a i _ 

RETURN 

e*#- 

c 

c 


« nmnnt  • 


IS 


CEUmi 

ccluips 

caiMi 

ccuiisi 

CCUISM 

ceutsti 

CEL1IS2I 

CCL1I93I 

ccliisii 

CELltSSI 

cattsM 

cmnn 

CEL1ISII 

ccliissi 
celumi 
cel  mu 

CEL1M2I 
CCtttlN 
CEillMI 
CELlfttSft 
CEtimi 
CELlfttfl 
CEL 111 II 
eeitun 
CELltMl 

ceumi 

CELlIf 21 
CftllMI 
CELllfM 

eetitrsi 

CEL2II1I 

.cEtmzi 

C CWMil 

celeiiii 

CEL2II9I 


COMPUTES  TKt  6CME4AU2G0  COMPLETE  ELLIPTIC  XNTC6RAL  OP 
"461111^4111^  ■ “ ■ 


CEL2IIM 


ft 

c 

- 

CALL 

CEL2(R£S,  AM,  A,ft,  IER) 

c 

DESCRIPT I ON  OP  PARAMETERS 

«*s 

c 

AM 

- MOOULUS  <I*?UTI 

c 

ft 

- FACTOR  OP  (JUAORATXC  T£RN  IN  \ 

c* 

I6MI  NO  ERROR 

c 

CEL2IIII 

CEtffttl 

CCL2II1I 


ss 


01 


c 

c 

c 
c 
c 
c 

c 

n c-  - 

c 

- ■*€“— 

__  C 

f-  c 

- - • c 

c_ 

c 

rr x 

c 

c 

• c 
I.  c 


FOR  AM  - ♦l.-l  TNI  RESULT  4kLJl  IS  SET  TO  1.ET9  IP  0 IS 
trtW  iP  • 10  MEtATUCi 


SPECIAL  CASiS  AH! 

MfrMP  PftTAINCft  trt*  t 

t«<l  OITAIMEO  JBTN  A • l,  I • CK •«  MMERE  CM  IS 
COMPtEMEMTAMT  MODULUS.  ~ “ 

• TMT  OWTAIMCU  WITH  A « I,  ft  ■ R 
OIMI  OITAIMEO  KITH  A a |,  ft  ■ l 

WHERE  M,  Tr'Of’  TT  OEPT ME  SPECIAL  CASES  Or  THE  6CNERILX2E3 
COMPLETE  ELLIPTIC  INTEGRAL  OP  ECONO  MXNO  IN  THE  USUAL 
~ NUTATION,  AND  IMI  HCUWLMT  T T TNC5E  POMCHOMS  MEINS 
TN£  MOOULUS. 

SUf ROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NOME  

METHOD — 

DEFINITION 

RES«INTlWRHT*t«*T*TlPrcaRTl  I I*T*T»»  IX*  <CM*T  »••!»>  • IX«T*T> 
SUMMED  OVER  T FROM  I TO  XNfXNXTT). 

rmorriDM - * 

LANOtNS  TRANSPORMATIOM  IS  USE)  FOR  CALCULATION. 

R.SULIRSCM,  * NUMERICAL  CALCULATION  OP  ELLIPTIC  XNTEtRAtS 
AMO  ILLIPTIC  FUNCTIONS-,  *(*3100*  SERIES  SPECIAL  FUNCTIONS, 
NUNLRISCHE  NATHLNATIM  VOL.  2,  1IIS,  PP.  21-SI. 


f - — 


SUIROUTINI  CtU«EfS,RMtR,l,ICRI 
lEMaft 

US  I MOOULUS 


6iO>l.-AM*IM 

“irmoir.m 

1 XtR*l 
"•RETURN 


Itt  WtUT  VALW  ■ OVERFLOW 


2 tPUIJiM 
I RtS>-X.E29 
'RETURN  • 

I RtS*l.U$ 


9 RiS*A 


CCUIUI 

CELEftllft 

CEL2I1SI 

C ft  Ml  It 

CEL2I1FI 

oetttit* 

CEL2I1SI 

CftWH 

CEL2I21I 

CCL2I2M 

CEL2I2SI 

cftfmi 

CEL2J2SI 

eatizii 

CEL2I2PI 

C£l2I2M 

CCL7I2SI 

CEL2IIII 

CEL2IJII 

CEL2I32I 

CEL2V3U 

CEL2I3II 

CEL2IITI 

CEL2I3M 

CEL2V32I 

CEL2I3I0 

CEL7TTTV 

CEL2IIII 

1CEL2IIU 

CEL2H2I 

CEL2II3I 

CEL2IIII 

ccimsi 

CEL2IIM 
CEL  71171 
CEL 2 I LI I 
CCUIISI 
•CEL2ISII 
CEL  21511 
CEL2I92I 
CEL 2 IS  3 1 
CEL2ISI0 
CCL2ISSI 
CEL20SI0 
CEL  21921 
CEL2ISII 
CEL 29991 
CEL2IMI 
CELZnil 
CEL2II2I 
CEL 2ft II 
CEL2IIII 
CEL2IISI 
CEL2tftftl 
CCL2II2I 
CELZtfttl 

cEimit 

CEL2I2U 
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2.  RcWN 

c 

C CUNfeJT*  I NUCHAL 

c 

* GtO»S9RTfft£OI 

2i  A*I«U 

AA«A 

IN*I*| 

U*l 

7 M«M*A A+G£0 
Jv  * — • •»*♦* 

AAsAN 

A«Rl*AKl 

A-<l*GtO«iftI 

Js  C 

- C-*  - US?  Of  AOOLMNC*  ~ 

c 

1HAAR1-GL0-I.t-«*AARIM,9,| 

• 6tO*SMT  <UcO#AARl> 

•*t  GtO*G£0»GtO 

GO  TO  7 

*-«rt*-.7t*JStT**A»r*Rt-  - ' ' ' 

Rf.  TURN 
INtJ 


CCL2I719 

CEL29729 

CCC29739 

CCL2I7AI 

ceinrM 
CEL29749 
CCL2I77I 
CCL29799 
C 0.297*9 
CEL29999 

cetmu 

CEL29929 
C 0.29939 
CCL2IINI 

cttmM 

cec2mi 

eetmn 

CEL29999 

CC12IM1 
cec 29111 

C2L29R19 

CEL29929 

tttmw 

CCL29*Nfi 

Cd.2t*99 


i 


r 
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